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

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)

eval:95:0: Type Checker: missing type for top-level

identifier;

 either undefined or missing a type annotation

  identifier: M16

  in: M

> (matrix-gram-schmidt M #t)

eval:96:0: Type Checker: missing type for top-level

identifier;

 either undefined or missing a type annotation

  identifier: M16

  in: #t

> (matrix-cols-orthogonal? (matrix-gram-schmidt M))

eval:97:0: Type Checker: missing type for top-level

identifier;

 either undefined or missing a type annotation

  identifier: M16

  in: M

> (matrix-orthonormal? (matrix-gram-schmidt M #t))

eval:98:0: Type Checker: missing type for top-level

identifier;

 either undefined or missing a type annotation

  identifier: M16

  in: #t

Examples with rank-deficient matrices:
> (matrix-gram-schmidt (matrix [[ 1 -2 1]
                                [-2  4 9]
                                [ 3 -6 7]]))

- : (Array Real)

(array #[#[1 5/7] #[-2 67/7] #[3 43/7]])

> (matrix-gram-schmidt (make-matrix 3 3 0))

- : (Array Real)

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

eval:102:0: Type Checker: missing type for top-level

identifier;

 either undefined or missing a type annotation

  identifier: M17

  in: #t

> (matrix-orthonormal?
   (matrix-gram-schmidt (matrix-gram-schmidt M) #t))

eval:103:0: Type Checker: missing type for top-level

identifier;

 either undefined or missing a type annotation

  identifier: M17

  in: #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))

eval:105:0: Type Checker: missing type for top-level

identifier;

 either undefined or missing a type annotation

  identifier: M18

  in: M

> (values Q R)

eval:106:0: Q19: unbound identifier;

 also, no #%top syntax transformer is bound

  in: Q19

> (matrix= M (matrix* Q R))

eval:107:0: Q19: unbound identifier;

 also, no #%top syntax transformer is bound

  in: Q19

The difference between full and reduced decompositions:
> (define M
    (matrix [[12 -51]
             [ 6 167]
             [-4  24]]))
> (define-values (Q1 R1) (matrix-qr M #f))

eval:109:0: Type Checker: missing type for top-level

identifier;

 either undefined or missing a type annotation

  identifier: M21

  in: #f

> (define-values (Q2 R2) (matrix-qr M #t))

eval:110:0: Type Checker: missing type for top-level

identifier;

 either undefined or missing a type annotation

  identifier: M21

  in: #t

> (values Q1 R1)

eval:111:0: Q122: unbound identifier;

 also, no #%top syntax transformer is bound

  in: Q122

> (values Q2 R2)

eval:112:0: Q224: unbound identifier;

 also, no #%top syntax transformer is bound

  in: Q224

> (matrix= M (matrix* Q1 R1))

eval:113:0: Q122: unbound identifier;

 also, no #%top syntax transformer is bound

  in: Q122

> (matrix= M (matrix* Q2 R2))

eval:114:0: Q224: unbound identifier;

 also, no #%top syntax transformer is bound

  in: Q224

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  ← prev  up  next →