Merge branch 'feat/sqrt-in-dml' of github.com:get4flo/systemds into feat/sqrt-in-dml

Author: trp-ex

hello.dml

@@ -0,0 +1,5 @@
+A = matrix(1, rows=2,cols=2)
+B = matrix(3, rows=2,cols=2)
+C = 10
+D = A %*% B + C * 2.1
+print( "D[1,1]:" + as.scalar(D[1,1]))

scripts/builtin/matrix_square_root.dml

@@ -0,0 +1,139 @@
+matrixSqrt = function(
+    Matrix[Double] X,
+    boolean iterMethod
+)return(
+    Matrix[Double] sqrt_x
+){
+    N = nrow(X);
+    D = ncol(X);
+
+    if( D == N ){
+        # Any non singualar square matrix has a square root
+        isDiag = isDiagonal(X)
+        if(isDiag) {
+            print("diag solution")
+            sqrt_x = sqrtDiagMatrix(X);
+        } else {
+            if(!iterMethod) {
+                #todo: check if all EigVal positive than possible
+                sqrt_x = sqrtEigInv(X);
+                print("eig inv solution")
+            } else {
+                # todo: iterative solution
+                #formular: (Denman–Beavers iteration)
+                print("iterative solution")
+                Y = X
+                #identity matrix
+                Z = diag(matrix(1.0, rows=N, cols=1))
+
+                for (x in 1:10) {
+                    Y_new = (1 / 2) * (Y + inv(Z))
+                    Z_new = (1 / 2) * (Z + inv(Y))
+                    Y = Y_new
+                    Z = Z_new
+                }
+                sqrt_x = Y
+            }
+        }
+    } else {
+        sqrt_x = matrix (0, rows=N, cols=D);
+    }
+}
+
+# assumes square and diagonal matrix
+sqrtDiagMatrix = function(
+        Matrix[Double] X
+)return(
+        Matrix[Double] sqrt_x
+){
+    N = nrow(X);
+
+    sqrt_x = matrix (0, rows=N, cols=N);
+    for (i in 1:N) {
+        value = X[i, i];
+        sqrt_x[i, i] = sqrt(value);
+    }
+}
+
+sqrtEigInv = function(
+        Matrix[Double] X
+)return(
+        Matrix[Double] sqrt_x
+){
+    [eValues, eVectors] = eigen(X);
+    # calculate X = VDV^(-1) -> S = sqrt(D) -> sqrt_x = VSV^(-1)
+    sqrtD = sqrtDiagMatrix(diag(eValues));
+    V_Inv = inv(eVectors);
+    sqrt_x = eVectors %*% sqrtD %*% V_Inv;
+}
+
+isDiagonal = function (
+    Matrix[Double] X
+)return(
+    boolean diagonal
+){
+    N = nrow(X);
+    D = ncol(X);
+    noCells = N * D;
+
+    diag = diag(diag(X));
+    compare = X == diag;
+    sameCells = sum(compare);
+
+    #all cells should be the same to be diagonal
+    diagonal = noCells == sameCells;
+}
+
+# testing area
+
+A = matrix(1, rows=2,cols=2)
+B = matrix(4, rows=2,cols=2)
+
+# easy test
+B[1,1] = 4.0
+B[1,2] = 0.0
+B[2,1] = 0.0
+B[2,2] = 4.0
+
+res = matrixSqrt(B, FALSE)
+print(as.scalar(res[1,1]))
+print(as.scalar(res[1,2]))
+print(as.scalar(res[2,1]))
+print(as.scalar(res[2,2]))
+
+# with this test for diag sqrt
+B[1,1] = 16.0
+B[1,2] = 21.0
+B[2,1] = 28.0
+B[2,2] = 37.0
+
+#matrixSqrt(A)
+#res = isDiagonal(A)
+#print(isDiagonal(A))
+res = matrixSqrt(B, FALSE)
+print(as.scalar(res[1,1]))
+print(as.scalar(res[1,2]))
+print(as.scalar(res[2,1]))
+print(as.scalar(res[2,2]))
+
+# iter test
+B[1,1] = 16.0
+B[1,2] = 21.0
+B[2,1] = 28.0
+B[2,2] = 37.0
+
+res = matrixSqrt(B, TRUE)
+print(as.scalar(res[1,1]))
+print(as.scalar(res[1,2]))
+print(as.scalar(res[2,1]))
+print(as.scalar(res[2,2]))
+
+#d = matrix("1.0 0.0
+#            0.0 1.0", 2, 2);
+#n = 3
+#I = diag(matrix(1.0, rows=n, cols=1))
+#print(toString(I, decimal=1))
+
+#test sqrt
+#C = sqrt(2)
+#print(C)