7.10 Row-Based Algorithms

procedure
(matrix-gauss-elim M
[ jordan?
unitize-pivot?
pivoting])
→ (Values (Matrix Number) (Listof Index))
  M : (Matrix Number)
  jordan? : Any = #f
  unitize-pivot? : Any = #f
  pivoting : (U 'first 'partial) = 'partial

Wikipedia: Gaussian elimination, Gauss-Jordan elimination Implements Gaussian elimination or Gauss-Jordan elimination.

If jordan? is true, row operations are done both above and below the pivot. If unitize-pivot? is true, the pivot’s row is scaled so that the pivot value is 1. When both are true, the algorithm is called Gauss-Jordan elimination, and the result matrix is in reduced row echelon form.

If pivoting is 'first, the first nonzero entry in the current column is used as the pivot. If pivoting is 'partial, the largest-magnitude nonzero entry is used, which improves numerical stability on average when M contains inexact entries.

The first return value is the result of Gaussian elimination.

The second return value is a list of indexes of columns that did not have a nonzero pivot.

See matrix-row-echelon for examples.

procedure
(matrix-row-echelon M
[ jordan?
unitize-pivot?
pivoting]) → (Matrix Number)
  M : (Matrix Number)
  jordan? : Any = #f
  unitize-pivot? : Any = #f
  pivoting : (U 'first 'partial) = 'partial

Wikipedia: Row echelon form Like matrix-gauss-elim, but returns only the result of Gaussian elimination.

Examples:

> (define M (matrix [[2 1 -1] [-3 -1 2] [-2 1 2]]))
> (matrix-row-echelon M)
- : #(struct:Array
      (Indexes Index (Boxof Boolean) (-> Void) (-> Indexes Real))
      #<syntax:.../array/typed-array-struct.rkt:56:13 prop:equal+hash>
      #<syntax:.../array/typed-array-struct.rkt:55:13 prop:custom-write>
      #<syntax:.../array/typed-array-struct.rkt:54:13 prop:custom-print-quotable>)
(mutable-array #[#[-3 -1 2] #[0 5/3 2/3] #[0 0 1/5]])
> (matrix-row-echelon M #t)
- : #(struct:Array
      (Indexes Index (Boxof Boolean) (-> Void) (-> Indexes Real))
      #<syntax:.../array/typed-array-struct.rkt:56:13 prop:equal+hash>
      #<syntax:.../array/typed-array-struct.rkt:55:13 prop:custom-write>
      #<syntax:.../array/typed-array-struct.rkt:54:13 prop:custom-print-quotable>)
(mutable-array #[#[-3 0 0] #[0 5/3 0] #[0 0 1/5]])
> (matrix-identity? (matrix-row-echelon M #t #t))
- : Boolean
#t

The last example shows that M is invertible.

Using matrix-row-echelon to solve a system of linear equations (without checking for failure):

> (define B (col-matrix [8 -11 -3]))
> (define MB (matrix-augment (list M B)))
> (matrix-col (matrix-row-echelon MB #t #t) 3)
- : #(struct:Array
      (Indexes Index (Boxof Boolean) (-> Void) (-> Indexes Real))
      #<syntax:.../array/typed-array-struct.rkt:56:13 prop:equal+hash>
      #<syntax:.../array/typed-array-struct.rkt:55:13 prop:custom-write>
      #<syntax:.../array/typed-array-struct.rkt:54:13 prop:custom-print-quotable>)
(array #[#[2] #[3] #[-1]])
> (matrix-solve M B)
- : #(struct:Array
      (Indexes Index (Boxof Boolean) (-> Void) (-> Indexes Real))
      #<syntax:.../array/typed-array-struct.rkt:56:13 prop:equal+hash>
      #<syntax:.../array/typed-array-struct.rkt:55:13 prop:custom-write>
      #<syntax:.../array/typed-array-struct.rkt:54:13 prop:custom-print-quotable>)
(array #[#[2] #[3] #[-1]])

Using matrix-row-echelon to invert a matrix (also without checking for failure):

> (define MI (matrix-augment (list M (identity-matrix 3))))
> (submatrix (matrix-row-echelon MI  #t #t) (::) (:: 3 #f))
- : #(struct:Array
      (Indexes Index (Boxof Boolean) (-> Void) (-> Indexes Real))
      #<syntax:.../array/typed-array-struct.rkt:56:13 prop:equal+hash>
      #<syntax:.../array/typed-array-struct.rkt:55:13 prop:custom-write>
      #<syntax:.../array/typed-array-struct.rkt:54:13 prop:custom-print-quotable>)
(array #[#[4 3 -1] #[-2 -2 1] #[5 4 -1]])
> (matrix-inverse M)
- : #(struct:Array
      (Indexes Index (Boxof Boolean) (-> Void) (-> Indexes Real))
      #<syntax:.../array/typed-array-struct.rkt:56:13 prop:equal+hash>
      #<syntax:.../array/typed-array-struct.rkt:55:13 prop:custom-write>
      #<syntax:.../array/typed-array-struct.rkt:54:13 prop:custom-print-quotable>)
(array #[#[4 3 -1] #[-2 -2 1] #[5 4 -1]])

procedure
(matrix-lu M [fail])
→ (Values (U F (Matrix Number)) (Matrix Number))
M : (Matrix Number)
fail : (-> F) = (λ () (error ...))

Wikipedia: LU decomposition Returns the LU decomposition of M (which must be a square-matrix?) if one exists. An LU decomposition exists if M can be put in row-echelon form without swapping rows.

Because matrix-lu returns a unit lower-triangular matrix (i.e. a lower-triangular matrix with only ones on the diagonal), the decomposition is unique if it exists.

Examples:

> (define-values (L U)
    (matrix-lu (matrix [[4 3] [6 3]])))
> (values L U)
- : (values
     #(struct:Array
       (Indexes Index (Boxof Boolean) (-> Void) (-> Indexes Real))
       #<syntax:.../array/typed-array-struct.rkt:56:13 prop:equal+hash>
       #<syntax:.../array/typed-array-struct.rkt:55:13 prop:custom-write>
       #<syntax:.../array/typed-array-struct.rkt:54:13 prop:custom-print-quotable>)
     #(struct:Array
       (Indexes Index (Boxof Boolean) (-> Void) (-> Indexes Real))
       #<syntax:.../array/typed-array-struct.rkt:56:13 prop:equal+hash>
       #<syntax:.../array/typed-array-struct.rkt:55:13 prop:custom-write>
       #<syntax:.../array/typed-array-struct.rkt:54:13 prop:custom-print-quotable>))
(mutable-array #[#[1 0] #[3/2 1]])
(mutable-array #[#[4 3] #[0 -3/2]])
> (matrix* L U)
- : #(struct:Array
      (Indexes Index (Boxof Boolean) (-> Void) (-> Indexes Real))
      #<syntax:.../array/typed-array-struct.rkt:56:13 prop:equal+hash>
      #<syntax:.../array/typed-array-struct.rkt:55:13 prop:custom-write>
      #<syntax:.../array/typed-array-struct.rkt:54:13 prop:custom-print-quotable>)
(array #[#[4 3] #[6 3]])

If M does not have an LU decomposition, the first result is the result of applying the failure thunk fail, and the second result is the original argument M:

> (matrix-lu (matrix [[0 1] [1 1]]))
matrix-lu: contract violation
  expected: LU-decomposable matrix
  given: (array #[#[0 1] #[1 1]])
> (matrix-lu (matrix [[0 1] [1 1]]) (λ () #f))
- : (values
     (U #(struct:Array
          (Indexes Index (Boxof Boolean) (-> Void) (-> Indexes Real))
          #<syntax:.../array/typed-array-struct.rkt:56:13 prop:equal+hash>
          #<syntax:.../array/typed-array-struct.rkt:55:13 prop:custom-write>
          #<syntax:.../array/typed-array-struct.rkt:54:13 prop:custom-print-quotable>)
        False)
     #(struct:Array
       (Indexes Index (Boxof Boolean) (-> Void) (-> Indexes Real))
       #<syntax:.../array/typed-array-struct.rkt:56:13 prop:equal+hash>
       #<syntax:.../array/typed-array-struct.rkt:55:13 prop:custom-write>
       #<syntax:.../array/typed-array-struct.rkt:54:13 prop:custom-print-quotable>))
#f
(array #[#[0 1] #[1 1]])

1	Constants and Elementary Functions
2	Flonums
3	Special Functions
4	Number Theory
5	Arbitrary-Precision Floating-Point Numbers (Bigfloats)
6	Arrays
7	Matrices and Linear Algebra
8	Statistics Functions
9	Probability Distributions
10	Stuff That Doesn’t Belong Anywhere Else

7.1	Introduction
7.2	Types, Predicates and Accessors
7.3	Construction
7.4	Conversion
7.5	Entrywise Operations and Arithmetic
7.6	Polymorphic Operations
7.7	Basic Operations
7.8	Inner Product Space Operations
7.9	Solving Systems of Equations
7.10	Row-Based Algorithms
7.11	Orthogonal Algorithms
7.12	Operator Norms and Comparing Matrices