incubator-systemml git commit: [SYSTEMML-1161] Recursive QR Factorization

[de...@apache.org]

Repository: incubator-systemml
Updated Branches:
  refs/heads/master 9cdce10d5 -> 4b6d468bb


[SYSTEMML-1161] Recursive QR Factorization

Implementation of QR Factorization using the Block Recursive algorithm
created by Elmroth and Gustavson. A simple representation of the algorithm
is presented by Golub (Algorithm 5.2.3).

Reference:
"Applying recursion to serial and parallel QR factorization leads to better
performance", E. Elmroth and F Gustavson, IBM J. RES. DEVELOP. VOL 44 JULY 2000
"Matrix Computations", Golub and Van Loan, 4th Edition, Chapter 5.

Closes #322.


Project: http://git-wip-us.apache.org/repos/asf/incubator-systemml/repo
Commit: http://git-wip-us.apache.org/repos/asf/incubator-systemml/commit/4b6d468b
Tree: http://git-wip-us.apache.org/repos/asf/incubator-systemml/tree/4b6d468b
Diff: http://git-wip-us.apache.org/repos/asf/incubator-systemml/diff/4b6d468b

Branch: refs/heads/master
Commit: 4b6d468bba4075ed0904d2818e1453bdf9b8a45e
Parents: 9cdce10
Author: Imran Younus <im...@gmail.com>
Authored: Sat Apr 29 09:23:25 2017 -0700
Committer: Deron Eriksson <de...@us.ibm.com>
Committed: Sat Apr 29 09:23:25 2017 -0700

----------------------------------------------------------------------
 scripts/staging/QR_recursive.dml | 90 +++++++++++++++++++++++++++++++++++
 1 file changed, 90 insertions(+)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/4b6d468b/scripts/staging/QR_recursive.dml
----------------------------------------------------------------------
diff --git a/scripts/staging/QR_recursive.dml b/scripts/staging/QR_recursive.dml
new file mode 100644
index 0000000..a3ee148
--- /dev/null
+++ b/scripts/staging/QR_recursive.dml
@@ -0,0 +1,90 @@
+#-------------------------------------------------------------
+#
+# Licensed to the Apache Software Foundation (ASF) under one
+# or more contributor license agreements.  See the NOTICE file
+# distributed with this work for additional information
+# regarding copyright ownership.  The ASF licenses this file
+# to you under the Apache License, Version 2.0 (the
+# "License"); you may not use this file except in compliance
+# with the License.  You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing,
+# software distributed under the License is distributed on an
+# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+# KIND, either express or implied.  See the License for the
+# specific language governing permissions and limitations
+# under the License.
+#
+#-------------------------------------------------------------
+
+# INPUT PARAMETERS:
+# -----------------------------------------------------------------------------
+# NAME   TYPE   DEFAULT  MEANING
+# -----------------------------------------------------------------------------
+# X      matrix  ---     input matrix to be factorized
+# k      Int     1000    smalles column block of input matrix
+# OUTDIR String  ---     output directory to save the results
+# OFMT   String  CSV     output format
+# -----------------------------------------------------------------------------
+
+# OUTPUT:
+# -----------------------------------------------------------------------------
+# NAME   TYPE    MEANING
+# -----------------------------------------------------------------------------
+# Q      matrix  Orthogonal matrix in thin QR factorization
+# R      matrix  Upper triangular matrix
+# -----------------------------------------------------------------------------
+
+X = read($X);
+k = ifdef($k, 1000);
+# k is used for the base case in the recursive algorithm.
+# If X has very large number of rows, choose k to be smaller
+# to fit the base case on single machine.
+ofmt = ifdef($OFMT, "CSV")
+
+QfromH = function(Matrix[double] H)
+  return(Matrix[double] Q) {
+    m = nrow(H);
+    n = ncol(H);
+    ones = matrix(1, m, 1);
+    eye = diag(ones);
+    Q = eye[,1:n];
+
+    for (j in n:1) {
+      v = H[j:m,j]
+      b = as.scalar(2/(t(v) %*% v))
+      Q[j:m, j:n] = Q[j:m, j:n] - (b * v) %*% (t(v) %*% Q[j:m, j:n])
+    }
+}
+
+QR_recursive = function(Matrix[double] A, int nb)
+  return(Matrix[double] Q, Matrix[double] R) {
+    n = ncol(A)
+
+    if (n <= nb) {
+      [H, R] = qr(A)
+      Q = QfromH(H)
+      R = R[1:n, 1:n]
+    }
+    else {
+      k = floor(n/2)
+      A1 = A[,1:k]
+      A2 = A[,k+1:n]
+
+      [Q1, R11] = QR_recursive(A1, nb)
+      R12 = t(Q1) %*% A2
+      A2 = A2 - Q1 %*% R12
+      [Q2, R22] = QR_recursive(A2, nb)
+
+      R21 = matrix(0, rows = nrow(R22), cols = ncol(R11))
+      Q = cbind(Q1, Q2)
+      R = rbind(cbind(R11, R12), cbind(R21, R22))
+    }
+}
+
+[Q, R] = QR_recursive(X, k)
+
+write(Q, $OUTDIR + "Q.csv", ofmt)
+write(R, $OUTDIR + "R.csv", ofmt)