Repository: systemml Updated Branches: refs/heads/master debe6d231 -> f8dc5c18b
[SYSTEMML-2270] Refactoring dml library for scalable decompositions This patch consolidates the individual dml files of the library for scalable decompositions into a single dml file which simplifies debugging and the inclusion in custom application scripts. Furthermore, the function Choleskey has been renamed to Cholesky. Project: http://git-wip-us.apache.org/repos/asf/systemml/repo Commit: http://git-wip-us.apache.org/repos/asf/systemml/commit/053e7ff0 Tree: http://git-wip-us.apache.org/repos/asf/systemml/tree/053e7ff0 Diff: http://git-wip-us.apache.org/repos/asf/systemml/diff/053e7ff0 Branch: refs/heads/master Commit: 053e7ff03369577dce88fcc82f93345faffb638d Parents: debe6d2 Author: Matthias Boehm <[email protected]> Authored: Sat Apr 21 19:50:24 2018 -0700 Committer: Matthias Boehm <[email protected]> Committed: Sat Apr 21 22:18:45 2018 -0700 ---------------------------------------------------------------------- scripts/staging/scalable_linalg/cholesky.dml | 45 ----- scripts/staging/scalable_linalg/inverse.dml | 35 ---- .../staging/scalable_linalg/linalg_decomp.dml | 186 +++++++++++++++++++ scripts/staging/scalable_linalg/lu.dml | 51 ----- scripts/staging/scalable_linalg/qr.dml | 61 ------ scripts/staging/scalable_linalg/solve.dml | 35 ---- .../staging/scalable_linalg/triangular_inv.dml | 80 -------- 7 files changed, 186 insertions(+), 307 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/systemml/blob/053e7ff0/scripts/staging/scalable_linalg/cholesky.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/cholesky.dml b/scripts/staging/scalable_linalg/cholesky.dml deleted file mode 100644 index 1bac1b2..0000000 --- a/scripts/staging/scalable_linalg/cholesky.dml +++ /dev/null @@ -1,45 +0,0 @@ -#------------------------------------------------------------- -# -# 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/systemml/blob/053e7ff0/scripts/staging/scalable_linalg/inverse.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/inverse.dml b/scripts/staging/scalable_linalg/inverse.dml deleted file mode 100644 index 112e351..0000000 --- a/scripts/staging/scalable_linalg/inverse.dml +++ /dev/null @@ -1,35 +0,0 @@ -#------------------------------------------------------------- -# -# 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/systemml/blob/053e7ff0/scripts/staging/scalable_linalg/linalg_decomp.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/linalg_decomp.dml b/scripts/staging/scalable_linalg/linalg_decomp.dml new file mode 100644 index 0000000..f665e34 --- /dev/null +++ b/scripts/staging/scalable_linalg/linalg_decomp.dml @@ -0,0 +1,186 @@ +#------------------------------------------------------------- +# +# 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. +# +#------------------------------------------------------------- + +Cholesky = 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 = Cholesky(A11, nb) + L11inv = U_triangular_inv(t(L11)) + L21 = A21 %*% L11inv + A22 = A22 - L21 %*% t(L21) + L22 = Cholesky(A22, nb) + L12 = matrix(0, nrow(L11), ncol(L22)) + + L = rbind(cbind(L11, L12), cbind(L21, L22)) + } +} + +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] = QR(A, nb) + Rinv = U_triangular_inv(R) + X = R %*% t(Q) +} + +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 = L_triangular_inv(L11) + U11inv = 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, nrow(A11), ncol(A22)) + Z21 = matrix(0, nrow(A22), 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)) + } +} + +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, nrow(R22), ncol(R11)) + + Q = cbind(Q1, Q2) + R = rbind(cbind(R11, R12), cbind(R21, R22)) + } +} + +# 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]) + } +} + +Solve = function(Matrix[double] A, Matrix[double] b, int nb) return(Matrix[double] x) { + # TODO: using thin QR may accumulate higher numerical errors. + # Some modifications as suggested by Golub may be explored + + [Q, R] = QR(A, nb) + Rinv = U_triangular_inv(R) + x = Rinv %*% t(Q) %*% b +} + +# 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, 2, 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, nrow(A11), 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, 2, 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, nrow(A22), ncol(A11)) + + A = rbind(cbind(A11, A12), cbind(A21, A22)) + } +} http://git-wip-us.apache.org/repos/asf/systemml/blob/053e7ff0/scripts/staging/scalable_linalg/lu.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/lu.dml b/scripts/staging/scalable_linalg/lu.dml deleted file mode 100644 index d1ad306..0000000 --- a/scripts/staging/scalable_linalg/lu.dml +++ /dev/null @@ -1,51 +0,0 @@ -#------------------------------------------------------------- -# -# 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/systemml/blob/053e7ff0/scripts/staging/scalable_linalg/qr.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/qr.dml b/scripts/staging/scalable_linalg/qr.dml deleted file mode 100644 index 7417a33..0000000 --- a/scripts/staging/scalable_linalg/qr.dml +++ /dev/null @@ -1,61 +0,0 @@ -#------------------------------------------------------------- -# -# 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/systemml/blob/053e7ff0/scripts/staging/scalable_linalg/solve.dml ---------------------------------------------------------------------- diff --git a/scripts/staging/scalable_linalg/solve.dml b/scripts/staging/scalable_linalg/solve.dml deleted file mode 100644 index 8dd07d2..0000000 --- a/scripts/staging/scalable_linalg/solve.dml +++ /dev/null @@ -1,35 +0,0 @@ -#------------------------------------------------------------- -# -# 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/systemml/blob/053e7ff0/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 deleted file mode 100644 index 75ff09e..0000000 --- a/scripts/staging/scalable_linalg/triangular_inv.dml +++ /dev/null @@ -1,80 +0,0 @@ -#------------------------------------------------------------- -# -# 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)) - } - }
