/** This is an incomplete class that implements matrix operations. Please feel free to implement missing functions! */ class Matrix { /** The array that will actually store the matrix data. Note that this is a raw 2-dimensional Frink array, and the indices are zero-based in contrast with most matrix notation which is one-based. The class should try to enforce that this is a rectangular 2-dimensional array, but the code doesn't currently enforce that with some constructors. */ var array /** The number of rows in the matrix. */ var rows /** The number of columns in the matrix. */ var cols /** Construct a new Matrix with the specified number of rows and columns and the specified default value for each element. */ new[rowCount, colCount, default=0] := { rows = rowCount cols = colCount array = makeArray[[rows, cols], default] } /** Construct a new Matrix from a rectangular 2-dimensional array. This currently does not verify that the array is rectangular. */ new[a is array] := { rows = length[a] cols = length[a@0] array = a } /** Sets the value of the specified element, using 1-based coordinates, as is common in matrix notation. Stupid mathematicians. */ set[row, col, val] := array@(row-1)@(col-1) = val /** Gets the value of the specified element, using 1-based coordinates, as is common in matrix notation. */ get[row, col] := array@(row-1)@(col-1) /** Multiply two matrices. The number of columns in column a should equal the number of rows in column b. The resulting matrix will have the same number of rows as matrix a and the same number of columns as matrix b. TODO: Use a faster algorithm for large arrays? This is O(n^3). TODO: Optimize this for small arrays by precalculating. */ class multiply[a, b] := { if a.cols != b.rows return undef else items = a.cols resRows = a.rows resCols = b.cols resultArray = makeArray[[resRows, resCols], 0] aa = a.array bb = b.array for row = 0 to resRows-1 for col = 0 to resCols-1 { sum = 0 for k = 0 to items-1 sum = sum + aa@row@k * bb@k@col resultArray@row@col = sum } return new Matrix[resultArray] } /** Multiply this matrix by the specified matrix and return the result. */ multiply[b] := multiply[this, b] /** Multiplies all elements of a matrix by a scalar. */ multiplyByScalar[s] := { a = makeArray[[rows, cols]] for rowNum = 0 to rows-1 { row = array@rowNum for colNum = 0 to cols-1 a@rowNum@colNum = row@colNum * s } return new Matrix[a] } /** Transposes the Matrix and returns the new transposed Matrix. This uses the built-in array.transpose[] method. */ transpose[] := new Matrix[array.transpose[]] /** Returns true if this Matrix is square, false otherwise. */ isSquare[] := rows == cols /** Exchange the specified-numbered rows in the matrix. The rows are 1-indexed to match standard matrix notation. */ exchangeRows[r1, r2] := { temp = array@(r2-1) array@(r2-1) = array@(r1-1) array@(r1-1) = temp } /** Exchange the specified-numbered columns in the matrix. The columns are 1-indexed to match standard matrix notation. */ exchangeCols[c1, c2] := { for r = 0 to rows-1 { temp = array@r@(c2-1) array@r@(c2-1) = array@r@(c1-1) array@r@(c1-1) = temp } } /** Returns a matrix as a string with rows separated with newlines. */ toString[] := { result = "" for r = 0 to rows-1 result = result + " " + array@r + "\n" return result } /** Returns the matrix formatted by formatTable with rows separated with newlines. */ toTable[] := { return formatTable[array] } /** Formats a matrix with Unicode box drawing characters. The result is a single string. */ formatMatrix[] := { return joinln[formatMatrixLines[]] } /** Returns the matrix formatted as a matrix with Unicode box drawing characters. The result is an array of lines that can be further formatted. */ formatMatrixLines[] := { m = formatTableLines[array] rows = length[m] width = length[m@0] result = new array result.push["\u250c" + repeat[" ", width] + "\u2510"] for n=0 to rows-1 result.push["\u2502" + m@n + "\u2502"] result.push["\u2514" + repeat[" ", width] + "\u2518"] return result } /** Gets the specified (1-based) column as a row array. */ getColumnAsArray[col] := { return deepCopy[array.getColumn[col-1]] } /** Gets the specified (1-based) row as a row array. */ getRowAsArray[row] := { return deepCopy[array@(row-1)] } /** Creates the LU (lower-upper) decomposition of the matrix. This creates two matrices, L and U, where L has the value 1 on the diagonal from top left to bottom right, like: 1 0 0 L = L21 1 0 L31 L32 1 U11 U12 U13 U = 0 U22 U23 0 0 U33 This is basically like solving equations by Gaussian elimination. */ LUDecompose[] := { return LUDecomposeCrout[] } /** This uses Crout's method to decompose a square matrix into an lower and upper triangular matrix. See: https://en.wikipedia.org/wiki/Crout_matrix_decomposition This creates two matrices, L and U, where U has the value 1 on the diagonal from top left to bottom right, like: L11 0 0 L = L21 L22 0 L31 L32 L33 1 U12 U13 U = 0 1 U23 0 0 1 returns [L, U] The original matrix can be obtained as L.multiply[U] Note that the Crout algorithm fails on certain matrices, for example m = new Matrix[[[1,1,1,-1],[1,1,-1,1],[1,-1,1,1],[-1,1,1,1]]] See: https://semath.info/src/inverse-cofactor-ex4.html for more on that matrix. https://en.wikipedia.org/wiki/Crout_matrix_decomposition */ LUDecomposeCrout[] := { n = rows L = new array[[rows, cols], 0] U = new array[[rows, cols], 0] sum = 0 for i=0 to n-1 U@i@i = 1 for j=0 to n-1 { for i=j to n-1 { sum = 0 for k=0 to j-1 sum = sum + L@i@k * U@k@j L@i@j = array@i@j - sum } for i=j to n-1 { sum = 0 for k=0 to j-1 sum = sum + L@j@k * U@k@i if L@j@j == 0 { println["Matrix.LUDecomposeCrout: det(L) is zero. Can't divide by zero."] return undef } U@j@i = (array@j@i - sum) / L@j@j } } return [new Matrix[L], new Matrix[U]] } /** Another algorithm for Crout LU decomposition. See: http://www.mosismath.com/Matrix_LU/Martix_LU.html This returns [L, U] as matrices which may be transposes of the matrices returned by LUDecomposeCrout[] . See its notes for more about this decomposition and its properties. */ LUDecomposeCrout2[] := { if ! isSquare[] { println["Matrix.LUDecomposeCrout2 only works on square arrays."] return undef } n = rows L = new array[[n,n], 0] U = new array[[n,n], 0] a = 0 for i = 0 to n-1 L@i@i = 1 for j = 0 to n-1 { for i = 0 to n-1 { if i >= j { U@j@i = array@j@i for k = 0 to j-1 { U@j@i = U@j@i - U@k@i * L@j@k } } if i > j { L@i@j = array@i@j for k = 0 to j-1 L@i@j = L@i@j - U@k@j * L@i@k diag = U@j@j if diag == 0 { println["Matrix.LUDecomposeCrout2: det(L) is zero. Can't divide by zero."] println[formatMatrix[L]] println[formatMatrix[U]] return undef } L@i@j = L@i@j / diag } } } return [new Matrix[L], new Matrix[U]] } /** This uses the Cholesky-Banachiewicz algorithm to decompose a square Hermitian matrix into a lower triangular matrix. If you transpose the lower triangular matrix L by L.transpose[], you get an upper triangular matrix which is symmetrical to the lower triangular matrix. Multiplying the lower triangular matrix by the upper triangular matrix, that is, L.multiply[L.transpose[]] gives you back the original matrix! See: https://en.m.wikipedia.org/wiki/Cholesky_decomposition */ CholeskyB[] := { if ! isHermitian[] { println["Matrix.CholeskyB only works on square Hermitian matrices."] return undef } n = rows L = new array[[n,n], 0] for i = 0 to n-1 for j = 0 to i { sum = 0 for k = 0 to j-1 sum = sum + L@i@k * L@j@k if i == j L@i@j = sqrt[array@i@i - sum] else L@i@j = (1 / L@j@j * (array@i@j -sum)) } return new Matrix[L] } /** This uses the Cholesky-Crout algorithm to decompose a square Hermitian matrix into a lower triangular matrix. If you transpose the lower triangular matrix L by L.transpose[], you get an upper triangular matrix which is symmetrical to the lower triangular matrix. Multiplying the lower triangular matrix by the upper triangular matrix, that is, L.multiply[L.transpose[]] gives you back the original matrix! See: https://en.m.wikipedia.org/wiki/Cholesky_decomposition */ CholeskyCrout[] := { if ! isHermitian[] { println["Matrix.CholeskyCrout only works on square Hermitian matrices."] return undef } n = rows L = new array[[n,n], 0] for j = 0 to n-1 { sum = 0 for k = 0 to j-1 sum = sum + (L@j@k)^2 L@j@j = sqrt[array@j@j - sum] for i = j+1 to n-1 { sum = 0 for k = 0 to j-1 sum = sum + L@i@k * L@j@k L@i@j = (1 / L@j@j * (array@i@j -sum)) } } return new Matrix[L] } /** This uses Doolittle's method to decompose a square matrix into an lower and upper triangular matrix. This creates two matrices, L and U, where L has the value 1 on the diagonal from top left to bottom right, like: 1 0 0 L = L21 1 0 L31 L32 1 U11 U12 U13 U = 0 U22 U23 0 0 U33 The original matrix can be obtained as L.multiply[U] Note that the Doolittle algorithm fails on certain matrices, for example m = new Matrix[[[1,1,1,-1],[1,1,-1,1],[1,-1,1,1],[-1,1,1,1]]] See: https://semath.info/src/inverse-cofactor-ex4.html for more on that matrix. See: https://www.geeksforgeeks.org/doolittle-algorithm-lu-decomposition/ */ LUDecomposeDoolittle[] := { if ! isSquare[] { println["Matrix.LUDecomposeDoolittle only works on square arrays."] return undef } n = rows L = new array[[n,n], 0] U = new array[[n,n], 0] // Decomposing matrix into Upper and Lower // triangular matrix for i = 0 to n-1 { // Upper Triangular for k = i to n-1 { // Summation of L(i, j) * U(j, k) sum = 0 for j = 0 to i-1 sum = sum + (L@i@j * U@j@k) // Evaluating U(i, k) U@i@k = array@i@k - sum } // Lower Triangular for k = i to n-1 { if i == k L@i@i = 1 // Diagonal as 1 else { // Summation of L(k, j) * U(j, i) sum = 0 for j = 0 to i-1 sum = sum + L@k@j * U@j@i diag = U@i@i if (diag == 0) { println["Matrix.LUDecomposeDoolittle: U@\$i@\$i is zero when solving for L@\$k@\$i. Can't divide by zero."] println[formatMatrix[L]] println[formatMatrix[U]] return undef } // Evaluating L(k, i) L@k@i = (array@k@i - sum) / diag } } } return [new Matrix[L], new Matrix[U]] } /** Returns the determinant of a square matrix. */ det[] := { if ! isSquare[] { println["Determinants are only defined for square matrices! Size was \$rows x \$cols"] return undef } // It's much faster to use precalculated determinant formulas for small // matrices than, to, say LUDecompose the things. if rows == 1 return array@0@0 if rows == 2 return array@0@0 array@1@1 - array@0@1 array@1@0 // ad - bc if rows == 3 return -array@0@2 array@1@1 array@2@0 + array@0@1 array@1@2 array@2@0 + array@0@2 array@1@0 array@2@1 - array@0@0 array@1@2 array@2@1 - array@0@1 array@1@0 array@2@2 + array@0@0 array@1@1 array@2@2 if rows == 4 return array@0@1 array@1@3 array@2@2 array@3@0 - array@0@1 array@1@2 array@2@3 array@3@0 - array@0@0 array@1@3 array@2@2 array@3@1 + array@0@0 array@1@2 array@2@3 array@3@1 - array@0@1 array@1@3 array@2@0 array@3@2 + array@0@0 array@1@3 array@2@1 array@3@2 + array@0@1 array@1@0 array@2@3 array@3@2 - array@0@0 array@1@1 array@2@3 array@3@2 + array@0@3 (array@1@2 array@2@1 array@3@0 - array@1@1 array@2@2 array@3@0 - array@1@2 array@2@0 array@3@1 + array@1@0 array@2@2 array@3@1 + array@1@1 array@2@0 array@3@2 - array@1@0 array@2@1 array@3@2) + (array@0@1 array@1@2 array@2@0 - array@0@0 array@1@2 array@2@1 - array@0@1 array@1@0 array@2@2 + array@0@0 array@1@1 array@2@2) array@3@3 + array@0@2 (-array@1@3 array@2@1 array@3@0 + array@1@1 array@2@3 array@3@0 + array@1@3 array@2@0 array@3@1 - array@1@0 array@2@3 array@3@1 - array@1@1 array@2@0 array@3@3 + array@1@0 array@2@1 array@3@3) if rows == 5 return array@0@2 array@1@4 array@2@3 array@3@1 array@4@0 - array@0@2 array@1@3 array@2@4 array@3@1 array@4@0 - array@0@1 array@1@4 array@2@3 array@3@2 array@4@0 + array@0@1 array@1@3 array@2@4 array@3@2 array@4@0 - array@0@2 array@1@4 array@2@1 array@3@3 array@4@0 + array@0@1 array@1@4 array@2@2 array@3@3 array@4@0 + array@0@2 array@1@1 array@2@4 array@3@3 array@4@0 - array@0@1 array@1@2 array@2@4 array@3@3 array@4@0 + array@0@2 array@1@3 array@2@1 array@3@4 array@4@0 - array@0@1 array@1@3 array@2@2 array@3@4 array@4@0 - array@0@2 array@1@1 array@2@3 array@3@4 array@4@0 + array@0@1 array@1@2 array@2@3 array@3@4 array@4@0 - array@0@2 array@1@4 array@2@3 array@3@0 array@4@1 + array@0@2 array@1@3 array@2@4 array@3@0 array@4@1 + array@0@0 array@1@4 array@2@3 array@3@2 array@4@1 - array@0@0 array@1@3 array@2@4 array@3@2 array@4@1 + array@0@2 array@1@4 array@2@0 array@3@3 array@4@1 - array@0@0 array@1@4 array@2@2 array@3@3 array@4@1 - array@0@2 array@1@0 array@2@4 array@3@3 array@4@1 + array@0@0 array@1@2 array@2@4 array@3@3 array@4@1 - array@0@2 array@1@3 array@2@0 array@3@4 array@4@1 + array@0@0 array@1@3 array@2@2 array@3@4 array@4@1 + array@0@2 array@1@0 array@2@3 array@3@4 array@4@1 - array@0@0 array@1@2 array@2@3 array@3@4 array@4@1 + array@0@1 array@1@4 array@2@3 array@3@0 array@4@2 - array@0@1 array@1@3 array@2@4 array@3@0 array@4@2 - array@0@0 array@1@4 array@2@3 array@3@1 array@4@2 + array@0@0 array@1@3 array@2@4 array@3@1 array@4@2 - array@0@1 array@1@4 array@2@0 array@3@3 array@4@2 + array@0@0 array@1@4 array@2@1 array@3@3 array@4@2 + array@0@1 array@1@0 array@2@4 array@3@3 array@4@2 - array@0@0 array@1@1 array@2@4 array@3@3 array@4@2 + array@0@1 array@1@3 array@2@0 array@3@4 array@4@2 - array@0@0 array@1@3 array@2@1 array@3@4 array@4@2 - array@0@1 array@1@0 array@2@3 array@3@4 array@4@2 + array@0@0 array@1@1 array@2@3 array@3@4 array@4@2 + array@0@2 array@1@4 array@2@1 array@3@0 array@4@3 - array@0@1 array@1@4 array@2@2 array@3@0 array@4@3 - array@0@2 array@1@1 array@2@4 array@3@0 array@4@3 + array@0@1 array@1@2 array@2@4 array@3@0 array@4@3 - array@0@2 array@1@4 array@2@0 array@3@1 array@4@3 + array@0@0 array@1@4 array@2@2 array@3@1 array@4@3 + array@0@2 array@1@0 array@2@4 array@3@1 array@4@3 - array@0@0 array@1@2 array@2@4 array@3@1 array@4@3 + array@0@1 array@1@4 array@2@0 array@3@2 array@4@3 - array@0@0 array@1@4 array@2@1 array@3@2 array@4@3 - array@0@1 array@1@0 array@2@4 array@3@2 array@4@3 + array@0@0 array@1@1 array@2@4 array@3@2 array@4@3 + array@0@2 array@1@1 array@2@0 array@3@4 array@4@3 - array@0@1 array@1@2 array@2@0 array@3@4 array@4@3 - array@0@2 array@1@0 array@2@1 array@3@4 array@4@3 + array@0@0 array@1@2 array@2@1 array@3@4 array@4@3 + array@0@1 array@1@0 array@2@2 array@3@4 array@4@3 - array@0@0 array@1@1 array@2@2 array@3@4 array@4@3 + array@0@4 (array@1@1 array@2@3 array@3@2 array@4@0 - array@1@1 array@2@2 array@3@3 array@4@0 - array@1@0 array@2@3 array@3@2 array@4@1 + array@1@0 array@2@2 array@3@3 array@4@1 - array@1@1 array@2@3 array@3@0 array@4@2 + array@1@0 array@2@3 array@3@1 array@4@2 + array@1@1 array@2@0 array@3@3 array@4@2 - array@1@0 array@2@1 array@3@3 array@4@2 + array@1@3 (array@2@2 array@3@1 array@4@0 - array@2@1 array@3@2 array@4@0 - array@2@2 array@3@0 array@4@1 + array@2@0 array@3@2 array@4@1 + array@2@1 array@3@0 array@4@2 - array@2@0 array@3@1 array@4@2) + (array@1@1 array@2@2 array@3@0 - array@1@0 array@2@2 array@3@1 - array@1@1 array@2@0 array@3@2 + array@1@0 array@2@1 array@3@2) array@4@3 + array@1@2 (-array@2@3 array@3@1 array@4@0 + array@2@1 array@3@3 array@4@0 + array@2@3 array@3@0 array@4@1 - array@2@0 array@3@3 array@4@1 - array@2@1 array@3@0 array@4@3 + array@2@0 array@3@1 array@4@3)) + (array@0@2 (-array@1@3 array@2@1 array@3@0 + array@1@1 array@2@3 array@3@0 + array@1@3 array@2@0 array@3@1 - array@1@0 array@2@3 array@3@1 - array@1@1 array@2@0 array@3@3 + array@1@0 array@2@1 array@3@3) + array@0@1 (array@1@3 array@2@2 array@3@0 - array@1@2 array@2@3 array@3@0 - array@1@3 array@2@0 array@3@2 + array@1@0 array@2@3 array@3@2 + array@1@2 array@2@0 array@3@3 - array@1@0 array@2@2 array@3@3) + array@0@0 (-array@1@3 array@2@2 array@3@1 + array@1@2 array@2@3 array@3@1 + array@1@3 array@2@1 array@3@2 - array@1@1 array@2@3 array@3@2 - array@1@2 array@2@1 array@3@3 + array@1@1 array@2@2 array@3@3)) array@4@4 + array@0@3 (-array@1@1 array@2@4 array@3@2 array@4@0 + array@1@1 array@2@2 array@3@4 array@4@0 + array@1@0 array@2@4 array@3@2 array@4@1 - array@1@0 array@2@2 array@3@4 array@4@1 + array@1@1 array@2@4 array@3@0 array@4@2 - array@1@0 array@2@4 array@3@1 array@4@2 - array@1@1 array@2@0 array@3@4 array@4@2 + array@1@0 array@2@1 array@3@4 array@4@2 + array@1@4 (-array@2@2 array@3@1 array@4@0 + array@2@1 array@3@2 array@4@0 + array@2@2 array@3@0 array@4@1 - array@2@0 array@3@2 array@4@1 - array@2@1 array@3@0 array@4@2 + array@2@0 array@3@1 array@4@2) - array@1@1 array@2@2 array@3@0 array@4@4 + array@1@0 array@2@2 array@3@1 array@4@4 + array@1@1 array@2@0 array@3@2 array@4@4 - array@1@0 array@2@1 array@3@2 array@4@4 + array@1@2 (array@2@4 array@3@1 array@4@0 - array@2@1 array@3@4 array@4@0 - array@2@4 array@3@0 array@4@1 + array@2@0 array@3@4 array@4@1 + array@2@1 array@3@0 array@4@4 - array@2@0 array@3@1 array@4@4)) // If we fell through to here, we have a larger matrix and have to try // to find its determinant through more inefficent methods. // TODO: Calculate determinant in terms of other determinants instead of // using LUDecompose when LUDecompose doesn't work. product = 1 [L, U] = LUDecomposeDoolittle[] // This will happen if the matrix is singular. if U == undef return undef // Multiply diagonals of the lower triangular matrix. The // upper triangular matrix has all ones on the diagonal. // Should there be a permutation matrix here to get the signs // right? for i=0 to rows-1 product = product * (U.array)@i@i return product } /** Create a square identity matrix with the specified number of rows and columns. The elements on the diagonal will be set to 1, the rest to zero. This requires Frink 2020-04-19 or later. */ class makeIdentity[dimension] := { return new Matrix[makeArray[[dimension, dimension], {|a,b| a==b ? 1 : 0}]] } /** Makes a diagonal matrix. This is passed in an array of elements (e.g. [1,2,3] that will make up the diagonal elements. This requires Frink 2020-04-22 or later. */ class makeDiagonal[array] := { d = length[array] return new Matrix[makeArray[[d,d], {|a,b,data| a==b ? data@a : 0}, array]] } /** Returns a matrix with the specified row and column removed. Row and column numbers are 1-indexed. This is used in many matrix operations including calculation of determinant or of the adjugate/adjoint matrix. */ removeRowColumn[rowToRemove, colToRemove] := { a = makeArray[[rows-1, cols-1]] newRow = 0 ROW: for origRow = 0 to rows-1 { if origRow == rowToRemove-1 next ROW newCol = 0 COL: for origCol = 0 to rows-1 { if origCol == colToRemove-1 next COL a@newRow@newCol = array@origRow@origCol newCol = newCol + 1 } newRow = newRow + 1 } return new Matrix[a] } /** Returns the adjugate or adjoint matrix. See: https://semath.info/src/inverse-cofactor-ex4.html TODO: Calculate hardcoded adjugate matrices for small matrices because this is expensive */ adjugate[] := { a = makeArray[[rows,cols]] for i = 0 to rows-1 for j = 0 to cols-1 a@i@j = (-1)^((i+1)+(j+1)) * removeRowColumn[j+1, i+1].det[] return new Matrix[a] } /** Returns the inverse of a matrix. See: https://semath.info/src/inverse-cofactor-ex4.html TODO: Calculate hardcoded inverse matrices for small matrices because this is expensive */ inverse[] := { if ! isSquare[] { println["Matrix.inverse only works on square arrays."] return undef } // It's much faster to calculate inverse matrices with hard-coded // equations for small matrices. if rows == 1 return 1 / array@0@0 if rows == 2 { invdet = 1/det[] // [ d -b ] // [ -c a ] / det return (new Matrix[[[array@1@1, -array@0@1], [-array@1@0, array@0@0]]]).multiplyByScalar[invdet] } if rows == 3 { invdet = 1/det[] return (new Matrix[[[-array@1@2 array@2@1 + array@1@1 array@2@2, array@0@2 array@2@1 - array@0@1 array@2@2, -array@0@2 array@1@1 + array@0@1 array@1@2], [ array@1@2 array@2@0 - array@1@0 array@2@2, -array@0@2 array@2@0 + array@0@0 array@2@2, array@0@2 array@1@0 - array@0@0 array@1@2], [-array@1@1 array@2@0 + array@1@0 array@2@1, array@0@1 array@2@0 - array@0@0 array@2@1, -array@0@1 array@1@0 + array@0@0 array@1@1]]]).multiplyByScalar[invdet] } if rows == 4 { invdet = 1/det[] return (new Matrix[[[-array@1@3 array@2@2 array@3@1 + array@1@2 array@2@3 array@3@1 + array@1@3 array@2@1 array@3@2 - array@1@1 array@2@3 array@3@2 - array@1@2 array@2@1 array@3@3 + array@1@1 array@2@2 array@3@3, array@0@3 array@2@2 array@3@1 - array@0@2 array@2@3 array@3@1 - array@0@3 array@2@1 array@3@2 + array@0@1 array@2@3 array@3@2 + array@0@2 array@2@1 array@3@3 - array@0@1 array@2@2 array@3@3, -array@0@3 array@1@2 array@3@1 + array@0@2 array@1@3 array@3@1 + array@0@3 array@1@1 array@3@2 - array@0@1 array@1@3 array@3@2 - array@0@2 array@1@1 array@3@3 + array@0@1 array@1@2 array@3@3, array@0@3 array@1@2 array@2@1 - array@0@2 array@1@3 array@2@1 - array@0@3 array@1@1 array@2@2 + array@0@1 array@1@3 array@2@2 + array@0@2 array@1@1 array@2@3 - array@0@1 array@1@2 array@2@3], [array@1@3 array@2@2 array@3@0 - array@1@2 array@2@3 array@3@0 - array@1@3 array@2@0 array@3@2 + array@1@0 array@2@3 array@3@2 + array@1@2 array@2@0 array@3@3 - array@1@0 array@2@2 array@3@3, -array@0@3 array@2@2 array@3@0 + array@0@2 array@2@3 array@3@0 + array@0@3 array@2@0 array@3@2 - array@0@0 array@2@3 array@3@2 - array@0@2 array@2@0 array@3@3 + array@0@0 array@2@2 array@3@3, array@0@3 array@1@2 array@3@0 - array@0@2 array@1@3 array@3@0 - array@0@3 array@1@0 array@3@2 + array@0@0 array@1@3 array@3@2 + array@0@2 array@1@0 array@3@3 - array@0@0 array@1@2 array@3@3, -array@0@3 array@1@2 array@2@0 + array@0@2 array@1@3 array@2@0 + array@0@3 array@1@0 array@2@2 - array@0@0 array@1@3 array@2@2 - array@0@2 array@1@0 array@2@3 + array@0@0 array@1@2 array@2@3], [-array@1@3 array@2@1 array@3@0 + array@1@1 array@2@3 array@3@0 + array@1@3 array@2@0 array@3@1 - array@1@0 array@2@3 array@3@1 - array@1@1 array@2@0 array@3@3 + array@1@0 array@2@1 array@3@3, array@0@3 array@2@1 array@3@0 - array@0@1 array@2@3 array@3@0 - array@0@3 array@2@0 array@3@1 + array@0@0 array@2@3 array@3@1 + array@0@1 array@2@0 array@3@3 - array@0@0 array@2@1 array@3@3, -array@0@3 array@1@1 array@3@0 + array@0@1 array@1@3 array@3@0 + array@0@3 array@1@0 array@3@1 - array@0@0 array@1@3 array@3@1 - array@0@1 array@1@0 array@3@3 + array@0@0 array@1@1 array@3@3, array@0@3 array@1@1 array@2@0 - array@0@1 array@1@3 array@2@0 - array@0@3 array@1@0 array@2@1 + array@0@0 array@1@3 array@2@1 + array@0@1 array@1@0 array@2@3 - array@0@0 array@1@1 array@2@3], [array@1@2 array@2@1 array@3@0 - array@1@1 array@2@2 array@3@0 - array@1@2 array@2@0 array@3@1 + array@1@0 array@2@2 array@3@1 + array@1@1 array@2@0 array@3@2 - array@1@0 array@2@1 array@3@2, -array@0@2 array@2@1 array@3@0 + array@0@1 array@2@2 array@3@0 + array@0@2 array@2@0 array@3@1 - array@0@0 array@2@2 array@3@1 - array@0@1 array@2@0 array@3@2 + array@0@0 array@2@1 array@3@2, array@0@2 array@1@1 array@3@0 - array@0@1 array@1@2 array@3@0 - array@0@2 array@1@0 array@3@1 + array@0@0 array@1@2 array@3@1 + array@0@1 array@1@0 array@3@2 - array@0@0 array@1@1 array@3@2, -array@0@2 array@1@1 array@2@0 + array@0@1 array@1@2 array@2@0 + array@0@2 array@1@0 array@2@1 - array@0@0 array@1@2 array@2@1 - array@0@1 array@1@0 array@2@2 + array@0@0 array@1@1 array@2@2]]]).multiplyByScalar[invdet] } // Otherwise use the adjugate function. return adjugate[].multiplyByScalar[1/det[]] } /** Returns the Kronecker product of a and b. https://en.wikipedia.org/wiki/Kronecker_product */ class KroneckerProduct[a is Matrix, b is Matrix] := { m = a.rows n = a.cols p = b.rows q = b.cols newRows = m p newCols = n q n = new array[[newRows, newCols], 0] for i=0 to newRows-1 for j=0 to newCols-1 n@i@j = (a.array)@(i div p)@(j div q) * (b.array)@(i mod p)@(j mod q) return new Matrix[n] } /** Return the Kronecker product of this matrix and b. */ KroneckerProduct[b is Matrix] := KroneckerProduct[this, b] /** Returns true if both matrices are equal (that is, they have the same dimensions and all array elements are equal.) */ equals[other is Matrix] := { return rows==other.rows and cols==other.cols and array==other.array } /** Returns a new Matrix where each element is the complex conjugate of the corresponding element in original Matrix. */ conjugate[] := { c = makeArray[[rows, cols], {|r,c,data| conjugate[data@r@c]}, array] return new Matrix[c] } /** Returns a new Matrix where the Matrix is first transposed and then each element is the complex conjugate of the corresponding element in the transposed Matrix. */ conjugateTranspose[] := { t = array.transpose[] c = makeArray[[rows, cols], {|r,c,data| conjugate[data@r@c]}, t] return new Matrix[c] } /** Returns true if this matrix is "Hermitian", or "self-adjoint", that is, the matrix must be square and equal to its own conjugate transpose. In other words, for all elements, array@i@j == conjugate[array@j@i] where conjugate is the complex conjugate of a number. */ isHermitian[] := { if rows != cols return false for i=0 to rows-1 for j = i+1 to rows-1 if array@i@j != conjugate[array@j@i] return false return true } /** Performs the QR decomposition of a matrix. This is based on the (real-only) implementation at: https://github.com/fiji/Jama/blob/master/src/main/java/Jama/QRDecomposition.java This is an internal implementation, primarily for use by the leastSquares method, that returns raw arrays [QR, rdiag]. If you're solving for least squares, use that routine directly instead of this. */ QRDecomposeInternal[] := { QR = deepCopy[array] // Copy the original matrix's array rdiag = new array[[cols], 0] // Stores the diagonal for k = 0 to cols-1 { // Compute 2-norm of k-th column nrm = 0 QR@0@k // Make units work right for i = k to rows-1 nrm = sqrt[nrm^2 + (QR@i@k)^2] if nrm != 0 QR@0@k // Make units work right { // Form k-th Householder vector if isNegativeUnit[QR@k@k] nrm = -nrm for i = k to rows-1 QR@i@k = QR@i@k / nrm QR@k@k = QR@k@k + 1 // Apply transformation to remaining columns for j = k+1 to cols-1 { s = undef for i=k to rows-1 { if s == undef s = QR@i@k * QR@i@j else s = s + QR@i@k * QR@i@j } s = -s / QR@k@k for i = k to rows-1 QR@i@j = QR@i@j + s * QR@i@k } } rdiag@k = -nrm } return [QR, rdiag] } /** This performs a QR decomposition of this matrix and returns [Q, R] as new matrices. Q is the orthogonal factor. R is the upper triangular factor. If you are doing a least-squares solve, call leastSquares directly instead. THINK ABOUT: Should this return the Householder vectors and rdiag? */ QRDecompose[] := { // The Q and R matrices are packed into the results of the following. [QR, rdiag] = QRDecomposeInternal[] // Extract the Q matrix Q = new array[[rows, cols], 0] for k = cols-1 to 0 step -1 { for i = 0 to rows-1 Q@i@k = 0 Q@k@k = 1 for j = k to cols-1 { if QR@k@k != 0 QR@k@k // Make units work right { s = 0 (QR@0@k) * Q@0@j // Make units work right for i = k to rows-1 s = s + QR@i@k * Q@i@j s = -s / QR@k@k for i = k to rows-1 Q@i@j = Q@i@j + s * QR@i@k } } } // Extract the R matrix R = new array[[cols, cols], 0] for i = 0 to cols-1 { for j = 0 to cols-1 { if i < j R@i@j = QR@i@j else if i == j R@i@j = rdiag@i else R@i@j = 0 } } return [new Matrix[Q], new Matrix[R]] } /** Return the least-squares solution (called X) of the system A * X = B where this Matrix is A and parameter B is B. This can be performed on an overdetermined system, where there are more measurements than equations. TODO: Try to solve exactly using x = A.inverse[] * b when the system is not overdetermined? See MatrixQRTest.frink for examples. This is the best discussion I've seen of least-squares fitting: https://www.aleksandrhovhannisyan.com/blog/the-method-of-least-squares/ https://www.aleksandrhovhannisyan.com/blog/least-squares-fitting/ */ leastSquares[B is Matrix] := { if rows != B.rows { println["Matrix.leastSquares: Matrix row dimensions must agree."] return undef } [QR, rdiag] = QRDecomposeInternal[] // Check if the matrix is "full rank," that is, that each row is // independent (not a multiple of) other rows. for j = 0 to cols-1 { if rdiag@j == 0 rdiag@j { println["Matrix.leastSquares: Row entries are not all independent."] return undef } } X = deepCopy[B.array] // Copy B's array // Compute Y = Q.transpose * B for k = 0 to cols-1 { for j = 0 to B.cols-1 { s = undef for i = k to rows-1 if s == undef s = QR@i@k * X@i@j else s = s + QR@i@k * X@i@j s = -s / QR@k@k for i = k to rows-1 X@i@j = X@i@j + s * QR@i@k } } // Solve R * X = y for k = cols-1 to 0 step -1 { for j = 0 to B.cols-1 X@k@j = X@k@j / rdiag@k for i = 0 to k-1 for j = 0 to B.cols-1 X@i@j = X@i@j - X@k@j * QR@i@k } X1 = new Matrix[X] return X1.getSubMatrix[0, cols-1, 0, B.cols-1] } /** Gets a submatrix of this matrix. */ getSubMatrix[fromRow, toRow, fromCol, toCol] := { a = new array[[toRow-fromRow+1, toCol-fromCol+1], 0] for row = fromRow to toRow for col = fromCol to toCol a@(row-fromRow)@(col-fromCol) = array@row@col // println["a is \$a"] return new Matrix[a] } /** Rounds values to the nearest integer if they are less than the specified relative error away from an integer and returns a new array. */ roundToInt[relerror = 1e-14] := { ret = new array[[rows, cols]] for r = 0 to rows-1 for c = 0 to cols-1 { v = array@r@c rnd = round[v] if abs[(rnd-v)/rnd] <= relerror ret@r@c = rnd else ret@r@c = v } return new Matrix[ret] } } "class Matrix loaded OK."