Download or view Matrix.frink in plain text format
/** 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
for i = k to rows-1
nrm = sqrt[nrm^2 + (QR@i@k)^2]
if nrm != 0
{
// 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 = 0
for i=k to rows-1
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
{
s = 0
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
{
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 = 0
for i = k to rows-1
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
}
// println["X is $X"]
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."
Download or view Matrix.frink in plain text format
This is a program written in the programming language Frink.
For more information, view the Frink
Documentation or see More Sample Frink Programs.
Alan Eliasen was born 19650 days, 7 hours, 42 minutes ago.