scripts/staging/scalable_linalg/linalg_decomp.dml

#-------------------------------------------------------------
#
# 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))
  }
}