7.11 Orthogonal Algorithms

procedure
(matrix-gram-schmidt M [normalize? start-col]) → (Array Number)
  M : (Matrix Number)
  normalize? : Any = #f
  start-col : Integer = 0

Wikipedia: Gram-Schmidt process Returns an array whose columns are orthogonal and span the same subspace as M’s columns. The number of columns in the result is the rank of M. If normalize? is true, the columns are also normalized.

Examples:

> (define M
    (matrix [[12 -51   4]
             [ 6 167 -68]
             [-4  24 -41]]))
> (matrix-gram-schmidt 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 #[#[12 -69 -58/5] #[6 158 6/5] #[-4 30 -33]])
> (matrix-gram-schmidt 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>)
(array #[#[6/7 -69/175 -58/175] #[3/7 158/175 6/175] #[-2/7 6/35 -33/35]])
> (matrix-cols-orthogonal? (matrix-gram-schmidt M))
- : Boolean
#t
> (matrix-orthonormal? (matrix-gram-schmidt M #t))
- : Boolean
#t

Examples with rank-deficient matrices:

> (matrix-gram-schmidt (matrix [[ 1 -2 1]
                                [-2  4 9]
                                [ 3 -6 7]]))
- : #(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 #[#[1 5/7] #[-2 67/7] #[3 43/7]])
> (matrix-gram-schmidt (make-matrix 3 3 0))
- : #(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 #[#[] #[] #[]])

When start-col is positive, the Gram-Schmidt process is begun on column start-col (but still using the previous columns to orthogonalize the remaining columns). This feature is generally not directly useful, but is used in the implementation of matrix-basis-extension.

* On the round-off error analysis of the Gram-Schmidt algorithm with reorthogonalization., Luc Giraud, Julien Langou and Miroslav Rozloznik. 2002. (PDF) While Gram-Schmidt with inexact matrices is known to be unstable, using it twice tends to remove instabilities:*

> (define M (matrix [[0.7 0.70711]
                     [0.70001 0.70711]]))
> (matrix-orthonormal?
   (matrix-gram-schmidt M #t))
- : Boolean
#f
> (matrix-orthonormal?
   (matrix-gram-schmidt (matrix-gram-schmidt M) #t))
- : Boolean
#t

procedure
(matrix-basis-extension M) → (Array Number)
M : (Matrix Number)

Returns additional orthogonal columns which, if augmented with M, would result in an orthogonal matrix of full rank. If M’s columns are normalized, the result’s columns are normalized.

procedure
(matrix-qr M) → (Values (Matrix Number) (Matrix Number))
  M : (Matrix Number)
(matrix-qr M full?) → (Values (Matrix Number) (Matrix Number))
  M : (Matrix Number)
  full? : Any

Wikipedia: QR decomposition Computes a QR-decomposition of the matrix M. The values returned are the matrices Q and R. If full? is #f, then a reduced decomposition is returned, otherwise a full decomposition is returned.

An orthonormal matrix has columns which are orthoginal, unit vectors. The (full) decomposition of a square matrix consists of two matrices: a orthogonal matrix Q and an upper triangular matrix R, such that QR = M.

For tall non-square matrices R, the triangular part of the full decomposition, contains zeros below the diagonal. The reduced decomposition leaves the zeros out. See the Wikipedia entry on QR decomposition for more details.

Examples:

> (define M
    (matrix [[12 -51   4]
             [ 6 167 -68]
             [-4  24 -41]]))
> (define-values (Q R) (matrix-qr M))
> (values Q R)
- : (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>))
(array #[#[6/7 -69/175 -58/175] #[3/7 158/175 6/175] #[-2/7 6/35 -33/35]])
(array #[#[14 21 -14] #[0 175 -70] #[0 0 35]])
> (matrix= M (matrix* Q R))
- : Boolean
#t

The difference between full and reduced decompositions:

> (define M
    (matrix [[12 -51]
             [ 6 167]
             [-4  24]]))
> (define-values (Q1 R1) (matrix-qr M #f))
> (define-values (Q2 R2) (matrix-qr M #t))
> (values Q1 R1)
- : (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>))
(array #[#[6/7 -69/175] #[3/7 158/175] #[-2/7 6/35]])
(array #[#[14 21] #[0 175]])
> (values Q2 R2)
- : (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>))
(array #[#[6/7 -69/175 58/175] #[3/7 158/175 -6/175] #[-2/7 6/35 33/35]])
(array #[#[14 21] #[0 175] #[0 0]])
> (matrix= M (matrix* Q1 R1))
- : Boolean
#t
> (matrix= M (matrix* Q2 R2))
- : Boolean
#t

The decomposition M = QR is useful for solving the equation Mx=v. Since the inverse of Q is simply the transpose of Q, Mx=v <=> QRx=v <=> Rx = Q^T v. And since R is upper triangular, the system can be solved by back substitution.

The algorithm used is Gram-Schmidt with reorthogonalization.

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