►Math Library
▼7 Matrices and Linear Algebra
On this page:
matrix-gram-schmidt
matrix-basis-extension
matrix-qr
  top  contents    up   

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 (Array Real) (Array Real))

(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 (Array Real) (Array Real))

(array #[#[6/7 -69/175] #[3/7 158/175] #[-2/7 6/35]])

(array #[#[14 21] #[0 175]])

> (values Q2 R2)

- : (values (Array Real) (Array Real))

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

  top  contents    up