Repository: incubator-systemml Updated Branches: refs/heads/master 5aa6fb75b -> 7ce130855
[SYSTEMML-1213] Scalable linear algebra DML implementations and accompanying DML tests for Cholesky, LU, QR, Solve, Inverse and util functions. Closes #368. Project: http://git-wip-us.apache.org/repos/asf/incubator-systemml/repo Commit: http://git-wip-us.apache.org/repos/asf/incubator-systemml/commit/7ce13085 Tree: http://git-wip-us.apache.org/repos/asf/incubator-systemml/tree/7ce13085 Diff: http://git-wip-us.apache.org/repos/asf/incubator-systemml/diff/7ce13085 Branch: refs/heads/master Commit: 7ce130855d66d94b79f1dabcb4eb584eccc9663c Parents: 5aa6fb7 Author: Imran Younus <imranyou...@gmail.com> Authored: Mon May 1 12:02:27 2017 -0700 Committer: Deron Eriksson <de...@us.ibm.com> Committed: Mon May 1 12:02:27 2017 -0700 ---------------------------------------------------------------------- scripts/staging/scalable_linalg/cholesky.dml | 45 +++++++++++ scripts/staging/scalable_linalg/inverse.dml | 35 +++++++++ scripts/staging/scalable_linalg/lu.dml | 51 +++++++++++++ scripts/staging/scalable_linalg/qr.dml | 61 +++++++++++++++ scripts/staging/scalable_linalg/solve.dml | 35 +++++++++ .../scalable_linalg/test/scalable_linalg | 1 + .../scalable_linalg/test/test_Cholesky.dml | 44 +++++++++++ .../staging/scalable_linalg/test/test_LU.dml | 44 +++++++++++ .../staging/scalable_linalg/test/test_QR.dml | 44 +++++++++++ .../staging/scalable_linalg/test/test_all.dml | 32 ++++++++ .../scalable_linalg/test/test_inverse.dml | 49 ++++++++++++ .../staging/scalable_linalg/test/test_solve.dml | 44 +++++++++++ .../test/test_triangular_inv.dml | 62 +++++++++++++++ scripts/staging/scalable_linalg/test/utils.dml | 38 ++++++++++ .../staging/scalable_linalg/triangular_inv.dml | 80 ++++++++++++++++++++ 15 files changed, 665 insertions(+) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/7ce13085/scripts/staging/scalable_linalg/cholesky.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/cholesky.dml b/scripts/staging/scalable_linalg/cholesky.dml new file mode 100644 index 0000000..1bac1b2 --- /dev/null +++ b/scripts/staging/scalable_linalg/cholesky.dml @@ -0,0 +1,45 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + +source("scalable_linalg/triangular_inv.dml") as inv + +Choleskey = function(Matrix[double] A, int nb) + return(Matrix[double] L) { + n = ncol(A) + + if (n <= nb) { + L = cholesky(A) + } else { + k = as.integer(floor(n/2)) + A11 = A[1:k,1:k] + A21 = A[k+1:n,1:k] + A22 = A[k+1:n,k+1:n] + + L11 = Choleskey(A11, nb) + L11inv = inv::U_triangular_inv(t(L11)) + L21 = A21 %*% L11inv + A22 = A22 - L21 %*% t(L21) + L22 = Choleskey(A22, nb) + L12 = matrix(0, rows=nrow(L11), cols=ncol(L22)) + + L = rbind(cbind(L11, L12), cbind(L21, L22)) + } +} http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/7ce13085/scripts/staging/scalable_linalg/inverse.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/inverse.dml b/scripts/staging/scalable_linalg/inverse.dml new file mode 100644 index 0000000..112e351 --- /dev/null +++ b/scripts/staging/scalable_linalg/inverse.dml @@ -0,0 +1,35 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + +source("scalable_linalg/triangular_inv.dml") as tinv +source("scalable_linalg/qr.dml") as decomp + +Inverse = function(Matrix[double] A, int nb) + return(Matrix[double] X) { + # TODO: Some recent papers discuss Block-Recursive algorithm for + # matrix inverse which can be explored instead of QR decomposition + + [Q, R] = decomp::QR(A, nb) + + Rinv = tinv::U_triangular_inv(R) + + X = R %*% t(Q) +} http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/7ce13085/scripts/staging/scalable_linalg/lu.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/lu.dml b/scripts/staging/scalable_linalg/lu.dml new file mode 100644 index 0000000..d1ad306 --- /dev/null +++ b/scripts/staging/scalable_linalg/lu.dml @@ -0,0 +1,51 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + +source("scalable_linalg/triangular_inv.dml") as inv + +LU = function(Matrix[double] A, int nb) + return(Matrix[double] P, Matrix[double] L, Matrix[double] U) { + n = ncol(A) + if (n <= nb) { + [P, L, U] = lu(A) + } else { + k = as.integer(floor(n/2)) + A11 = A[1:k,1:k] + A12 = A[1:k,k+1:n] + A21 = A[k+1:n,1:k] + A22 = A[k+1:n,k+1:n] + + [P11, L11, U11] = LU(A11, nb) + L11inv = inv::L_triangular_inv(L11) + U11inv = inv::U_triangular_inv(U11) + U12 = L11inv %*% (P11 %*% A12) + A22 = A22 - A21 %*% U11inv %*% U12 + [P22, L22, U22] = LU(A22, nb) + L21 = P22 %*% A21 %*% U11inv + + Z12 = matrix(0, rows=nrow(A11), cols=ncol(A22)) + Z21 = matrix(0, rows=nrow(A22), cols=ncol(A11)) + + L = rbind(cbind(L11, Z12), cbind(L21, L22)) + U = rbind(cbind(U11, U12), cbind(Z21, U22)) + P = rbind(cbind(P11, Z12), cbind(Z21, P22)) + } +} http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/7ce13085/scripts/staging/scalable_linalg/qr.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/qr.dml b/scripts/staging/scalable_linalg/qr.dml new file mode 100644 index 0000000..7417a33 --- /dev/null +++ b/scripts/staging/scalable_linalg/qr.dml @@ -0,0 +1,61 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + +# calculate Q from Housholder matrix +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 = 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(A1, nb) + R12 = t(Q1) %*% A2 + A2 = A2 - Q1 %*% R12 + [Q2, R22] = QR(A2, nb) + R21 = matrix(0, rows = nrow(R22), cols = ncol(R11)) + + Q = cbind(Q1, Q2) + R = rbind(cbind(R11, R12), cbind(R21, R22)) + } +} http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/7ce13085/scripts/staging/scalable_linalg/solve.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/solve.dml b/scripts/staging/scalable_linalg/solve.dml new file mode 100644 index 0000000..8dd07d2 --- /dev/null +++ b/scripts/staging/scalable_linalg/solve.dml @@ -0,0 +1,35 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + +source("scalable_linalg/triangular_inv.dml") as inv +source("scalable_linalg/qr.dml") as decomp + +Solve = function(Matrix[double] A, Matrix[double] b, int nb) + return(Matrix[double] x) { + # TODO: using thin QR may accumulate higher numerical erros. + # Some modifications as suggested by Golub may be explored + + [Q, R] = decomp::QR(A, nb) + + Rinv = inv::U_triangular_inv(R) + + x = Rinv %*% t(Q) %*% b +} http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/7ce13085/scripts/staging/scalable_linalg/test/scalable_linalg ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/test/scalable_linalg b/scripts/staging/scalable_linalg/test/scalable_linalg new file mode 120000 index 0000000..b0e1ab9 --- /dev/null +++ b/scripts/staging/scalable_linalg/test/scalable_linalg @@ -0,0 +1 @@ +../../scalable_linalg \ No newline at end of file http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/7ce13085/scripts/staging/scalable_linalg/test/test_Cholesky.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/test/test_Cholesky.dml b/scripts/staging/scalable_linalg/test/test_Cholesky.dml new file mode 100644 index 0000000..f602e93 --- /dev/null +++ b/scripts/staging/scalable_linalg/test/test_Cholesky.dml @@ -0,0 +1,44 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + +source("scalable_linalg/cholesky.dml") as decomp + +test_Cholesky = function() { + print("Testing Choleskey Decomposition.") + n = 1000 + b = 100 # smallest block size + eps = n*n*1e-12 # for lareger matrices eps should be larger + + # create a symmetric random matrix + A = rand(rows=n, cols=n, min=-1.0, max=1.0, pdf="uniform") + X = t(A) %*% A + + L = decomp::Choleskey(X, b) + + # check if X = LL^T. Infinity norm of (X - LL^T) must be close to zero + diff = X - L %*% t(L) + sup_norm = max(abs(diff)) + if (sup_norm > eps) { + print("ERROR: Cholesky decomposition does not reproduce original matrix") + } +} + +tmp = test_Cholesky() http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/7ce13085/scripts/staging/scalable_linalg/test/test_LU.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/test/test_LU.dml b/scripts/staging/scalable_linalg/test/test_LU.dml new file mode 100644 index 0000000..307535e --- /dev/null +++ b/scripts/staging/scalable_linalg/test/test_LU.dml @@ -0,0 +1,44 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + +source("scalable_linalg/lu.dml") as decomp + +test_LU = function() { + print("Testing LU Decomposition.") + n = 1000 + b = 100 # smallest block size + eps = n*n*1e-12 # for lareger matrices eps should be larger + + # create random square matrix + A = rand(rows=n, cols=n, min=-1.0, max=1.0, pdf="uniform") + + [P, L, U] = decomp::LU(A, b) + + # check if PA = LU. Infinity norm of (PA - LU) must be close to zero + diff = P %*% A - L %*% U + sup_norm = max(abs(diff)) + print(sup_norm) + if (sup_norm > eps) { + print("ERROR: LU decomposition does not reproduce original matrix") + } +} + +tmp = test_LU() http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/7ce13085/scripts/staging/scalable_linalg/test/test_QR.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/test/test_QR.dml b/scripts/staging/scalable_linalg/test/test_QR.dml new file mode 100644 index 0000000..83f2214 --- /dev/null +++ b/scripts/staging/scalable_linalg/test/test_QR.dml @@ -0,0 +1,44 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + +source("scalable_linalg/qr.dml") as decomp + +test_QR = function() { + print("Testing QR Decomposition.") + n = 1000 + m = 2000 + b = 100 # smallest block size + eps = n*m*1e-12 # for lareger matrices eps should be larger + + # create a symmetric random matrix + A = rand(rows=m, cols=n, min=-1.0, max=1.0, pdf="uniform") + + [Q, R] = decomp::QR(A, b) + + # check if A = QR. Infinity norm of (A - QR) must be close to zero + diff = A - Q %*% R + sup_norm = max(abs(diff)) + if (sup_norm > eps) { + print("ERROR: QR decomposition does not reproduce original matrix") + } +} + +tmp = test_QR() http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/7ce13085/scripts/staging/scalable_linalg/test/test_all.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/test/test_all.dml b/scripts/staging/scalable_linalg/test/test_all.dml new file mode 100644 index 0000000..bba391c --- /dev/null +++ b/scripts/staging/scalable_linalg/test/test_all.dml @@ -0,0 +1,32 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + +source("test_inverse.dml") as test1 +source("test_LU.dml") as test2 +source("test_QR.dml") as test3 +source("test_Cholesky.dml") as test4 + +print("Testing Scalable Linalg") + +tmp = test1::test_inverse() +tmp = test2::test_LU() +#tmp = test3::test_QR() +tmp = test4::test_Cholesky() http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/7ce13085/scripts/staging/scalable_linalg/test/test_inverse.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/test/test_inverse.dml b/scripts/staging/scalable_linalg/test/test_inverse.dml new file mode 100644 index 0000000..2bd0f5d --- /dev/null +++ b/scripts/staging/scalable_linalg/test/test_inverse.dml @@ -0,0 +1,49 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + +source("scalable_linalg/inverse.dml") as inv + +test_inverse = function() { + print("Testing inverse of square matrices") + n = 1000 + b = 100 # smallest block size + eps = n*n*1e-12 + + A = rand(rows=n, cols=n, min=-1.0, max=1.0, pdf="uniform") + + Ainv = inv::Inverse(A, b) + + # this product should produce identity matrix + AAinv = A %*% Ainv + + # create identity matrix + ones = matrix(1.0, n, 1); + I = diag(ones); + + # check if AA^-1 = I. Infinity norm of (I - AA^-1) must be close to zero + diff = I - AAinv + sup_norm = max(abs(diff)) + if (sup_norm > eps) { + print("ERROR: inverse of square matrix fails") + } +} + +tmp = test_inverse() http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/7ce13085/scripts/staging/scalable_linalg/test/test_solve.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/test/test_solve.dml b/scripts/staging/scalable_linalg/test/test_solve.dml new file mode 100644 index 0000000..f741c0f --- /dev/null +++ b/scripts/staging/scalable_linalg/test/test_solve.dml @@ -0,0 +1,44 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + +source("scalable_linalg/solve.dml") as lsq + +test_solve = function() { + print("Testing linear Solve function.") + n = 1000 + m = 2000 + b = 100 + eps = 1e-9 + + # generate data for test + A = rand(rows=m, cols=n, min=-2, max=2) + x = rand(rows=n, cols=1, min=-5, max=5) + y = A %*% x + x1 = lsq::Solve(A, y, b) + + diff = abs(x - x1) + sup_norm = max(diff) + if (sup_norm > eps) { + print("ERROR: least squares solve fails") + } +} + +tmp = test_solve() http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/7ce13085/scripts/staging/scalable_linalg/test/test_triangular_inv.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/test/test_triangular_inv.dml b/scripts/staging/scalable_linalg/test/test_triangular_inv.dml new file mode 100644 index 0000000..13d730a --- /dev/null +++ b/scripts/staging/scalable_linalg/test/test_triangular_inv.dml @@ -0,0 +1,62 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + +source("scalable_linalg/triangular_inv.dml") as inv + +test_inverse = function() { + print("Testing inverse of triangular matrices") + n = 1000 + eps = n*n*1e-12 + + A = rand(rows=n, cols=n, min=-1.0, max=1.0, pdf="uniform") + + # create lower/upper triangular matrices using lu deomposition + [P,L,U] = lu(A) + + # calculate inverse + Linv = inv::L_triangular_inv(L) + Uinv = inv::U_triangular_inv(U) + + # these products should produce identity matrices + LLinv = L %*% Linv + UUinv = U %*% Uinv + + # create identity matrix + ones = matrix(1.0, n, 1); + I = diag(ones); + + # check if LL^-1 = I. Infinity norm of (I - LL^-1) must be close to zero + diff = I - LLinv + sup_norm = max(abs(diff)) + if (sup_norm > eps) { + print("ERROR: inverse of lower triangular matrix fails") + } + + # check if UU^-1 = I. Inifinity norm of (I - UU^-1) must be close to zero + diff = I - UUinv + sup_norm = max(abs(diff)) + + if (sup_norm > eps) { + print("ERROR: inverse of upper triangular matrix fails") + } +} + +tmp = test_inverse() http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/7ce13085/scripts/staging/scalable_linalg/test/utils.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/test/utils.dml b/scripts/staging/scalable_linalg/test/utils.dml new file mode 100644 index 0000000..765c97a --- /dev/null +++ b/scripts/staging/scalable_linalg/test/utils.dml @@ -0,0 +1,38 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + + +rel_error = function(double true_val, double approx) + return(double err) { + err = abs(true_val - approx)/max(abs(true_val), 1e-9) + } + +check_equal = function(double x1, double x2, double eps) + return (boolean eq) { + eq = TRUE + diff = abs(x1 - x2) + largest = max(abs(x2), abs(x2)) + + if (diff > largest*eps) { + print("ERROR: vlaues not equal: " + x1 + " != " + x2 ) + eq = FALSE + } + } http://git-wip-us.apache.org/repos/asf/incubator-systemml/blob/7ce13085/scripts/staging/scalable_linalg/triangular_inv.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/triangular_inv.dml b/scripts/staging/scalable_linalg/triangular_inv.dml new file mode 100644 index 0000000..75ff09e --- /dev/null +++ b/scripts/staging/scalable_linalg/triangular_inv.dml @@ -0,0 +1,80 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + + +# Inverse of lower triangular matrix +L_triangular_inv = function(Matrix[double] L) + return(Matrix[double] A) { + n = ncol(L) + + if (n == 1) { + A = 1/L[1,1] + } else if (n == 2) { + A = matrix(0, rows=2, cols=2) + A[1,1] = L[2,2] + A[2,2] = L[1,1] + A[2,1] = -L[2,1] + A = A/(as.scalar(L[1,1] * L[2,2])) + } else { + k = as.integer(floor(n/2)) + + L11 = L[1:k,1:k] + L21 = L[k+1:n,1:k] + L22 = L[k+1:n,k+1:n] + + A11 = L_triangular_inv(L11) + A22 = L_triangular_inv(L22) + A12 = matrix(0, rows=nrow(A11), cols=ncol(A22)) + A21 = -A22 %*% L21 %*% A11 + + A = rbind(cbind(A11, A12), cbind(A21, A22)) + } + } + +# Inverse of upper triangular matrix +U_triangular_inv = function(Matrix[double] U) + return(Matrix[double] A) { + n = ncol(U) + + if (n == 1) { + A = 1/U[1,1] + } else if (n == 2) { + A = matrix(0, rows=2, cols=2) + A[1,1] = U[2,2] + A[2,2] = U[1,1] + + A[1,2] = -U[1,2] + A = A/(as.scalar(U[1,1] * U[2,2])) + } else { + k = as.integer(floor(n/2)) + + U11 = U[1:k,1:k] + U12 = U[1:k,k+1:n] + U22 = U[k+1:n,k+1:n] + + A11 = U_triangular_inv(U11) + A22 = U_triangular_inv(U22) + A12 = -A11 %*% U12 %*% A22 + A21 = matrix(0, rows=nrow(A22), cols=ncol(A11)) + + A = rbind(cbind(A11, A12), cbind(A21, A22)) + } + }
