Matrix.frink

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 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."


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 20217 days, 15 hours, 57 minutes ago.