Finished matrix documentation, attendant fixes; please merge to 5.3.2

* Narrowed type of `submatrix' to only sensible argument types

* `matrix-invertible?' now returns #f when given a non-square matrix
  instead of raising an error

* Allowed `matrix-diagonal' to operate on non-square matrices
(cherry picked from commit 24561e25e4)
This commit is contained in:
Neil Toronto 2013-01-17 19:06:51 -07:00 committed by Ryan Culpepper
parent 8dd7ed2004
commit 6708c385e6
5 changed files with 492 additions and 168 deletions

View File

@ -70,7 +70,8 @@
[else
(unsafe-array-ref a ((inst vector Index) i j))]))
(: submatrix (All (A) (Matrix A) Slice-Spec Slice-Spec -> (Matrix A)))
(: submatrix (All (A) (Matrix A) (U Slice (Sequenceof Integer)) (U Slice (Sequenceof Integer))
-> (Array A)))
(define (submatrix a row-range col-range)
(array-slice-ref (ensure-matrix 'submatrix a) (list row-range col-range)))
@ -120,11 +121,11 @@
(: matrix-diagonal (All (A) ((Matrix A) -> (Array A))))
(define (matrix-diagonal a)
(define m (square-matrix-size a))
(define-values (m n) (matrix-shape a))
(define proc (unsafe-array-proc a))
(array-default-strict
(unsafe-build-array
((inst vector Index) m)
((inst vector Index) (fxmin m n))
(λ: ([js : Indexes])
(define i (unsafe-vector-ref js 0))
(proc ((inst vector Index) i i))))))
@ -305,8 +306,11 @@
(: matrix-trace (case-> ((Matrix Real) -> Real)
((Matrix Number) -> Number)))
(define (matrix-trace a)
(parameterize ([array-strictness #f])
(array-all-sum (matrix-diagonal a))))
(cond [(square-matrix? a)
(parameterize ([array-strictness #f])
(array-all-sum (matrix-diagonal a)))]
[else
(raise-argument-error 'matrix-trace "square-matrix?" a)]))
;; ===================================================================================================
;; Row/column operations

View File

@ -69,7 +69,8 @@
(: matrix-invertible? ((Matrix Number) -> Boolean))
(define (matrix-invertible? M)
(not (zero? (matrix-determinant M))))
(and (square-matrix? M)
(not (zero? (matrix-determinant M)))))
(: matrix-inverse (All (A) (case-> ((Matrix Real) -> (Matrix Real))
((Matrix Real) (-> A) -> (U A (Matrix Real)))

View File

@ -6,7 +6,7 @@
racket/list
math plot
(only-in typed/racket/base
ann inst : λ: define: make-predicate ->
ann inst : λ: define: make-predicate -> case-> Nothing
Flonum Real Boolean Any Integer Index Natural Exact-Positive-Integer
Nonnegative-Real Sequenceof Fixnum Values Number Float-Complex
All U List Vector Listof Vectorof Struct FlVector
@ -41,9 +41,16 @@ For now, if you need speed, use the @racketmodname[typed/racket] language.
@defmodule[math/matrix]
Documentation for this module is currently under construction.
Intro topics: definitions, case-> types, non-strictness
Like all of @racketmodname[math], @racketmodname[math/matrix] is a work in progress.
Most of the basic algorithms are implemented, but some are still in planning.
Possibly the most useful unimplemented algorithms are
@itemlist[@item{LUP decomposition (currently, LU decomposition is implemented, in @racket[matrix-lu])}
@item{@racket[matrix-solve] for upper-triangular matrices}
@item{Singular value decomposition (SVD)}
@item{Eigendecomposition}
@item{Decomposition-based solvers}
@item{Pseudoinverse, least-squares fitting}]
All of these are planned for the next Racket release.
@local-table-of-contents[]
@ -51,6 +58,76 @@ Intro topics: definitions, case-> types, non-strictness
@;{==================================================================================================}
@section[#:tag "matrix:intro"]{Introduction}
From the point of view of the functions in @racketmodname[math/matrix], a @deftech{matrix} is an
@racket[Array] with two axes and at least one entry, or an array for which @racket[matrix?] returns
@racket[#t].
Technically, a matrix's entries may be any type, and some fully polymorphic matrix functions such as
@racket[matrix-row] and @racket[matrix-map] operate on any kind of matrix.
Other functions, such as @racket[matrix+], require their input matrices to contain either
@racket[Real] or @racket[Number] values.
@subsection[#:tag "matrix:function-types"]{Function Types}
The documentation for @racketmodname[math/matrix] functions use the type @racket[Matrix], a synonym
of @racket[Array], when the function either requires that an argument is a @tech{matrix} or ensures
that a return value is a matrix.
Most functions that implement matrix algorithms are documented as accepting @racket[(Matrix Number)]
values. This includes @racket[(Matrix Real)], which is a subtype. Most of these functions have a more
precise type than is documented. For example, @racket[matrix-conjugate] actually has the type
@racketblock[(case-> ((Matrix Real) -> (Matrix Real))
((Matrix Number) -> (Matrix Number)))]
even though it is documented as having the less precise type
@racket[((Matrix Number) -> (Matrix Number))].
Precise function types allow Typed Racket to prove more facts about @racketmodname[math/matrix]
client programs. In particular, it is usually easy for it to prove that matrix expressions on real
matrices return real matrices:
@interaction[#:eval typed-eval
(matrix-conjugate (matrix [[1 2 3] [4 5 6]]))]
@subsection[#:tag "matrix:failure"]{Failure Arguments}
In many matrix operations, such as inversion, failure is easy to detect during computation, but
is just as expensive to detect ahead of time as the operation itself.
In these cases, the functions implementing the operations accept an optional @deftech{failure thunk},
or a zero-argument function that returns the result of the operation in case of failure.
For example, the (simplified) type of @racket[matrix-inverse] is
@racketblock[(All (F) (case-> ((Matrix Number) -> (Matrix Number))
((Matrix Number) (-> F) -> (U F (Matrix Number)))))]
Thus, if a failure thunk is given, the call site is required to check for return values of type
@racket[F] explicitly.
The default failure thunk, which raises an error, has type @racket[(-> Nothing)].
For such failure thunks, @racket[(U F (Matrix Number))] is equivalent to @racket[(Matrix Number)],
because @racket[Nothing] is part of every type. (In Racket, any expression may raise an error.)
Thus, in this case, no explicit test for values of type @racket[F] is necessary (though of course they
may be caught using @racket[with-handlers] or similar).
@subsection[#:tag "matrix:broadcasting"]{Broadcasting}
Pointwise matrix operations do not @tech{broadcast} their arguments when given matrices with
different sizes:
@interaction[#:eval typed-eval
(matrix+ (identity-matrix 2) (matrix [[10]]))]
If you need broadcasting, use array operations:
@interaction[#:eval typed-eval
(array+ (identity-matrix 2) (matrix [[10]]))]
@subsection{Strictness}
Functions exported by @racketmodname[math/matrix] return @tech{strict} or @tech{nonstrict} arrays
based on the value of the @racket[array-strictness] parameter.
See @secref{array:nonstrict} for details.
@;{==================================================================================================}
@section[#:tag "matrix:types"]{Types, Predicates and Accessors}
@defform[(Matrix A)]{
@ -59,7 +136,7 @@ Equivalent to @racket[(Array A)], but used for values @racket[M] for which @rack
}
@defproc[(matrix? [arr (Array A)]) Boolean]{
Returns @racket[#t] when @racket[arr] is a @deftech{matrix}: a nonempty array with exactly two axes.
Returns @racket[#t] when @racket[arr] is a @tech{matrix}: a nonempty array with exactly two axes.
@examples[#:eval typed-eval
(matrix? (array 10))
(matrix? (array #[1 2 3]))
@ -161,7 +238,9 @@ Returns a matrix with @racket[xs] along the diagonal and @racket[0] everywhere e
The length of @racket[xs] must be positive.
}
@margin-note{@hyperlink["http://en.wikipedia.org/wiki/Block_matrix#Block_diagonal_matrices"]{Wikipedia: Block-diagonal matrices}}
@define[block-diagonal-url]{http://en.wikipedia.org/wiki/Block_matrix#Block_diagonal_matrices}
@margin-note*{@hyperlink[block-diagonal-url]{Wikipedia: Block-diagonal matrices}}
@defproc[(block-diagonal-matrix [Xs (Listof (Matrix A))]) (Matrix (U A 0))]{
Returns a matrix with matrices @racket[Xs] along the diagonal and @racket[0] everywhere else.
The length of @racket[Xs] must be positive.
@ -172,7 +251,9 @@ The length of @racket[Xs] must be positive.
(row-matrix [4 5 6])))]
}
@margin-note{@hyperlink["http://en.wikipedia.org/wiki/Vandermonde_matrix"]{Wikipedia: Vandermonde matrix}}
@define[vandermonde-url]{http://en.wikipedia.org/wiki/Vandermonde_matrix}
@margin-note*{@hyperlink[vandermonde-url]{Wikipedia: Vandermonde matrix}}
@defproc[(vandermonde-matrix [xs (Listof Number)] [n Integer]) (Matrix Number)]{
Returns an @racket[m]×@racket[n] Vandermonde matrix, where @racket[m = (length xs)].
@examples[#:eval typed-eval
@ -313,11 +394,11 @@ where each entry in @racket[M] is multiplied with @racket[z].
(matrix-scale (matrix [[1 2] [3 4]]) 2)]
}
@defproc*[([(matrix-map [f (A -> R)] [arr0 (Matrix A)]) (Matrix R)]
[(matrix-map [f (A B Ts ... -> R)] [arr0 (Matrix A)] [arr1 (Matrix B)] [arrs (Matrix Ts)]
@defproc*[([(matrix-map [f (A -> R)] [M (Matrix A)]) (Matrix R)]
[(matrix-map [f (A B Ts ... -> R)] [M0 (Matrix A)] [M1 (Matrix B)] [N (Matrix Ts)]
...)
(Matrix R)])]{
Like @racket[array-map], but requires at least one array argument and never @tech{broadcasts}.
Like @racket[array-map], but requires at least one matrix argument and never @tech{broadcasts}.
@examples[#:eval untyped-eval
(matrix-map sqr (matrix [[1 2] [3 4]]))
(matrix-map + (matrix [[1 2] [3 4]])
@ -344,47 +425,52 @@ tolerant to floating-point error.
@defproc[(matrix-ref [M (Matrix A)] [i Integer] [j Integer]) A]{
Returns the entry on row @racket[i] and column @racket[j].
@examples[#:eval untyped-eval
(define A (matrix ([1 2 3] [4 5 6])))
(matrix-ref A 0 2)
(matrix-ref A 1 2)]
}
@defthing[submatrix Procedure]{
@;{ TODO
(: submatrix (All (A) (Matrix A) Slice-Spec Slice-Spec -> (Matrix A)))
(define (submatrix a row-range col-range)
(array-slice-ref (ensure-matrix 'submatrix a) (list row-range col-range)))
}
TODO
(define M (matrix ([1 2 3] [4 5 6])))
(matrix-ref M 0 2)
(matrix-ref M 1 2)]
}
@deftogether[(@defproc[(matrix-row [M (Matrix A)] [i Integer]) (Matrix A)]
@defproc[(matrix-col [M (Matrix A)] [j Integer]) (Matrix A)])]{
Returns the given row or column.
Returns the @racket[i]th row or @racket[j]th column as a matrix.
@examples[#:eval untyped-eval
(define A (matrix ([1 2 3] [4 5 6])))
(matrix-row A 1)
(matrix-col A 0)]
(define M (matrix ([1 2 3] [4 5 6])))
(matrix-row M 1)
(matrix-col M 0)]
}
@defproc[(submatrix [M (Matrix A)]
[is (U Slice (Sequenceof Integer))]
[js (U Slice (Sequenceof Integer))])
(Array A)]{
Returns a submatrix or subarray of @racket[M], where @racket[is] and @racket[js] specify
respectively the rows and columns to keep. Like @racket[array-slice-ref], but constrained so the
result has exactly two axes.
@examples[#:eval typed-eval
(submatrix (identity-matrix 5) (:: 1 #f 2) (::))
(submatrix (identity-matrix 5) '() '(1 2 4))]
Note that @racket[submatrix] may return an empty array, which is not a @tech{matrix}.
}
@defproc[(matrix-diagonal [M (Matrix A)]) (Array A)]{
Returns array of the elements on the diagonal of the square matrix.
Returns array of the entries on the diagonal of @racket[M].
@examples[#:eval untyped-eval
(define A (matrix ([1 2 3] [4 5 6] [7 8 9])))
(matrix-diagonal A)]
(matrix-diagonal
(matrix ([1 2 3] [4 5 6] [7 8 9])))]
}
@deftogether[(@defproc[(matrix-upper-triangle [M (Matrix A)]) (Matrix A)]
@defproc[(matrix-lower-triangle [M (Matrix A)]) (Matrix A)])]{
@deftogether[(@defproc[(matrix-upper-triangle [M (Matrix A)]) (Matrix (U A 0))]
@defproc[(matrix-lower-triangle [M (Matrix A)]) (Matrix (U A 0))])]{
The function @racket[matrix-upper-triangle] returns an upper
triangular matrix (entries below the diagonal are zero) with
elements from the given matrix. Likewise the function
@racket[matrix-lower-triangle] returns an lower triangular
entries from the given matrix. Likewise the function
@racket[matrix-lower-triangle] returns a lower triangular
matrix.
@examples[#:eval untyped-eval
(define A (matrix ([1 2 3] [4 5 6] [7 8 9])))
(matrix-upper-triangle A)
(matrix-lower-triangle A)]
(define M (array+ (array 1) (axis-index-array #(5 7) 1)))
M
(matrix-upper-triangle M)
(matrix-lower-triangle M)]
}
@deftogether[(@defproc[(matrix-rows [M (Matrix A)]) (Listof (Matrix A))]
@ -392,45 +478,61 @@ matrix.
The functions respectively returns a list of the rows or columns
of the matrix.
@examples[#:eval untyped-eval
(define A (matrix ([1 2 3] [4 5 6])))
(matrix-rows A)
(matrix-cols A)]
(define M (matrix ([1 2 3] [4 5 6])))
(matrix-rows M)
(matrix-cols M)]
}
@deftogether[(@defproc[(matrix-augment [Ms (Listof (Matrix A))]) (Matrix A)]
@defproc[(matrix-stack [Ms (Listof (Matrix A))]) (Matrix A)])]{
The function @racket[matrix-augment] returns a matrix whose columns are
the columns of the matrices in @racket[Ms]. This implies that the matrices
the columns of the matrices in @racket[Ms]. The matrices
in list must have the same number of rows.
The function @racket[matrix-stack] returns a matrix whose rows are
the rows of the matrices in @racket[Ms]. This implies that the matrices
the rows of the matrices in @racket[Ms]. The matrices
in list must have the same number of columns.
@examples[#:eval untyped-eval
(define A (matrix ([1 1] [1 1])))
(define B (matrix ([2 2] [2 2])))
(define C (matrix ([3 3] [3 3])))
(matrix-augment (list A B C))
(matrix-stack (list A B C))]
(define M0 (matrix ([1 1] [1 1])))
(define M1 (matrix ([2 2] [2 2])))
(define M2 (matrix ([3 3] [3 3])))
(matrix-augment (list M0 M1 M2))
(matrix-stack (list M0 M1 M2))]
}
@deftogether[
(@defproc[(matrix-map-rows
[f ((Matrix A) -> (U #f (Matrix B)))] [M (Matrix A)] [fail (-> F) (λ () #f)]) (Matrix B)]
@defproc[(matrix-map-cols
[f ((Matrix A) -> (U #f (Matrix B)))] [M (Matrix A)] [fail (-> F) (λ () #f)]) (Matrix B)])]{
In the simple case the function @racket[matrix-map-rows] applies the function @racket[f]
to each row of @racket[M]. If the rows are called @racket[r0], @racket[r1], ... then
the result matrix has the rows @racket[(f r0)], @racket[(f r1)], ... .
In the three argument case, the result of @racket[(fail)] is used,
if @racket[f] returns @racket[#f].
The function @racket[matrix-map-cols] works likewise but on rows.
(@defproc*[([(matrix-map-rows
[f ((Matrix A) -> (Matrix B))] [M (Matrix A)]) (Matrix B)]
[(matrix-map-rows
[f ((Matrix A) -> (U #f (Matrix B)))] [M (Matrix A)] [fail (-> F)])
(U F (Matrix B))])])]{
The two-argument case applies the function @racket[f] to each row of @racket[M].
If the rows are called @racket[r0], @racket[r1], ..., the result matrix has the rows
@racket[(f r0)], @racket[(f r1)], ....
@examples[#:eval untyped-eval
(define A (matrix ([1 2 3] [4 5 6] [7 8 9] [10 11 12])))
(define M (matrix ([1 2 3] [4 5 6] [7 8 9] [10 11 12])))
(define (double-row r) (matrix-scale r 2))
(matrix-map-rows double-row A)]
(matrix-map-rows double-row M)]
In the three argument case, if @racket[f] returns @racket[#f], the result of @racket[(fail)] is
returned:
@interaction[#:eval typed-eval
(define Z (make-matrix 4 4 0))
Z
(matrix-map-rows (λ: ([r : (Matrix Real)])
(matrix-normalize r 2 (λ () #f)))
Z
(λ () 'FAILURE))]
}
@deftogether[
(@defproc*[([(matrix-map-cols
[f ((Matrix A) -> (Matrix B))] [M (Matrix A)]) (Matrix B)]
[(matrix-map-cols
[f ((Matrix A) -> (U #f (Matrix B)))] [M (Matrix A)] [fail (-> F)])
(U F (Matrix B))])])]{
Like @racket[matrix-map-cols], but maps @racket[f] over columns.
}
@ -439,16 +541,16 @@ The function @racket[matrix-map-cols] works likewise but on rows.
@section[#:tag "matrix:basic"]{Basic Operations}
@defproc[(matrix-conjugate [M (Matrix A)]) (Matrix A)]{
Returns a matrix where each element of the given matrix is conjugated.
@defproc[(matrix-conjugate [M (Matrix Number)]) (Matrix Number)]{
Returns a matrix where each entry of the given matrix is conjugated.
@examples[#:eval untyped-eval
(matrix-conjugate (matrix ([1 +i] [-1 2+i])))]
}
@margin-note{@hyperlink["http://en.wikipedia.org/wiki/Transpose"]{Wikipedia: Transpose}}
@margin-note*{@hyperlink["http://en.wikipedia.org/wiki/Transpose"]{Wikipedia: Transpose}}
@deftogether[(@defproc[(matrix-transpose [M (Matrix A)]) (Matrix A)]
@defproc[(matrix-hermitian [M (Matrix A)]) (Matrix A)])]{
@margin-note{@hyperlink["http://en.wikipedia.org/wiki/Hermitian_matrix"]{Wikipedia: Hermitian}}
@defproc[(matrix-hermitian [M (Matrix Number)]) (Matrix Number)])]{
@margin-note*{@hyperlink["http://en.wikipedia.org/wiki/Hermitian_matrix"]{Wikipedia: Hermitian}}
Returns the transpose or the hermitian of the matrix.
The hermitian of a matrix is the conjugate of the transposed matrix.
For a real matrix these operations return the the same result.
@ -457,71 +559,68 @@ For a real matrix these operations return the the same result.
(matrix-hermitian (matrix ([1 +i] [2 +2i] [3 +3i])))]
}
@margin-note{@hyperlink["http://en.wikipedia.org/wiki/Trace_(linear_algebra)"]{Wikipedia: Trace}}
@defproc[(matrix-trace [M (Matrix Number)]) (Matrix Number)]{
@margin-note*{@hyperlink["http://en.wikipedia.org/wiki/Trace_(linear_algebra)"]{Wikipedia: Trace}}
@defproc[(matrix-trace [M (Matrix Number)]) Number]{
Returns the trace of the square matrix. The trace of matrix is the
the sum of the diagonal elements.
the sum of the diagonal entries.
@examples[#:eval untyped-eval
(matrix-trace (matrix ([1 2] [3 4])))]
}
@;{==================================================================================================}
@section[#:tag "matrix:inner"]{Inner Product Space Operations}
The set of matrices of a given size forms a vector space.
Since the vector space of @racket[mx1] matrices is isomorphic to
the vector space of vectors of size @racket[m], any inner
product (or norm) induce an inner product (or norm) on vectors.
Put differently, the following innner products and norm, even
though defined on general matrices, work on vectors in the
form of column and row matrices.
The following functions treat matrices as vectors in an inner product space.
It often makes most sense to use these vector-space functions only for
@tech[#:key "row matrix"]{row matrices} and @tech[#:key "column matrix"]{column matrices}, which are
essentially vectors as we normally think of them.
There are exceptions, however, such as the fact that the Frobenius or Euclidean norm (implemented by
@racket[matrix-2norm]) can be used to measure error between matrices in a way that meets certain
reasonable criteria (specifically, it is submultiplicative).
See @secref{matrix:op-norm} for similar functions (e.g. norms and angles) defined by considering
matrices as operators between inner product spaces consisting of column matrices.
@margin-note{@hyperlink["http://en.wikipedia.org/wiki/Norm_(mathematics)"]{Wikipedia: Norm}}
@deftogether[(@defproc[(matrix-1norm [M (Matrix Number)]) Number]
@defproc[(matrix-2norm [M (Matrix Number)]) Number]
@defproc[(matrix-inf-norm [M (Matrix Number)]) Number]
@defproc*[([(matrix-norm [M (Matrix Number)]) Number]
[(matrix-norm [M (Matrix Number)] [p Real]) Number])])]{
The first three functions compute the L1-norm, the L2-norm, and, the L∞-norm respectively.
@margin-note*{@hyperlink["http://en.wikipedia.org/wiki/Norm_(mathematics)"]{Wikipedia: Norm}}
@deftogether[(@defproc[(matrix-1norm [M (Matrix Number)]) Nonnegative-Real]
@defproc[(matrix-2norm [M (Matrix Number)]) Nonnegative-Real]
@defproc[(matrix-inf-norm [M (Matrix Number)]) Nonnegative-Real]
@defproc[(matrix-norm [M (Matrix Number)] [p Real 2]) Nonnegative-Real])]{
Respectively compute the L@subscript{1} norm, L@subscript{2} norm, L@subscript{∞}, and
L@subscript{p} norm.
The L1-norm is also known under the names Manhattan- or taxicab-norm.
The L1-norm of a matrix is the sum of magnitudes of the entries in the matrix.
The L@subscript{1} norm is also known under the names Manhattan or taxicab norm.
The L@subscript{1} norm of a matrix is the sum of magnitudes of the entries in the matrix.
The L2-norm is also known under the names Euclidean- or Frobenius-norm.
The L2-norm of a matrix is the square root of the sum of squares of
The L@subscript{2} norm is also known under the names Euclidean or Frobenius norm.
The L@subscript{2} norm of a matrix is the square root of the sum of squares of
magnitudes of the entries in the matrix.
The L∞-norm is also known as the maximum- or infinity-norm.
The L∞-norm computes the maximum magnitude of the entries in the matrix.
The L@subscript{∞} norm is also known as the maximum or infinity norm.
The L@subscript{∞} norm computes the maximum magnitude of the entries in the matrix.
The function @racket[matrix-norm] computes the Lp-norm.
For a number @racket[p>=1] the @racket[p]th root of the sum
of all entries to the @racket[p]th power.
For @racket[p >= 1], @racket[matrix-norm] computes the L@subscript{@racket[p]} norm:
the @racket[p]th root of the sum of all entry magnitudes to the @racket[p]th power.
@;{MathJax would be nice to have in Scribble...}
If no @racket[p] is given, the 2-norm (Eucledian) is used.
@examples[#:eval untyped-eval
(matrix-1norm (col-matrix [1 2]))
(matrix-2norm (col-matrix [1 2]))
(matrix-inf-norm (col-matrix [1 2]))
(matrix-norm (col-matrix [1 2]) 3)]
(matrix-norm (col-matrix [1 2]) 3)
(matrix-norm (col-matrix [1 2]) +inf.0)]
}
@defproc*[([(matrix-dot [M (Matrix Number)]) Nonnegative-Real]
[(matrix-dot [M1 (Matrix Number)] [M2 (Matrix Number)]) Number])]{
[(matrix-dot [M (Matrix Number)] [N (Matrix Number)]) Number])]{
The call @racket[(matrix-dot M1 M2)] computes the Frobenius inner product of the
The call @racket[(matrix-dot M N)] computes the Frobenius inner product of the
two matrices with the same shape.
In other words the sum of @racket[(* a (conjugate b))] is computed where
@racket[a] runs over the entries in @racket[M1] and @racket[b] runs over
the corresponding entries in @racket[M2].
@racket[a] runs over the entries in @racket[M] and @racket[b] runs over
the corresponding entries in @racket[N].
The call @racket[(matrix-dot M)] computes @racket[(matrix-dot M M)] efficiently.
@examples[#:eval untyped-eval
@ -529,79 +628,130 @@ The call @racket[(matrix-dot M)] computes @racket[(matrix-dot M M)] efficiently.
(+ (* 1 3) (* 2 4))]
}
@defproc[(matrix-cos-angle [M (Matrix Number)]) Number]{
@defproc[(matrix-cos-angle [M (Matrix Number)] [N (Matrix Number)]) Number]{
Returns the cosine of the angle between two matrices w.r.t. the inner produce space induced by
the Frobenius inner product. That is it returns
@racket[(/ (matrix-dot M N) (* (matrix-2norm M) (matrix-2norm N)))]
@examples[#:eval untyped-eval
(define e1 (col-matrix [1 0]))
(define e2 (col-matrix [0 1]))
(matrix-cos-angle e1 e2)
(matrix-cos-angle e1 (matrix+ e1 e2))]
(define M (col-matrix [1 0]))
(define N (col-matrix [0 1]))
(matrix-cos-angle M N)
(matrix-cos-angle M (matrix+ M N))]
}
@defproc[(matrix-angle [M1 (Matrix Number)] [M2 (Matrix Number)]) Number]{
Equivalent to @racket[(acos (matrix-cos-angle M0 M1))].
@defproc[(matrix-angle [M (Matrix Number)] [N (Matrix Number)]) Number]{
Equivalent to @racket[(acos (matrix-cos-angle M N))].
@examples[#:eval untyped-eval
(require racket/math) ; for radians->degrees
(define e1 (col-matrix [1 0]))
(define e2 (col-matrix [0 1]))
(radians->degrees (matrix-angle e1 e2))
(radians->degrees (matrix-angle e1 (matrix+ e1 e2)))]
(require (only-in math/base radians->degrees))
(define M (col-matrix [1 0]))
(define N (col-matrix [0 1]))
(radians->degrees (matrix-angle M N))
(radians->degrees (matrix-angle M (matrix+ M N)))]
}
@defproc[(matrix-normalize [M (Matrix Number)] [p Real 2] [fail (-> A) raise-argument-error]) (U (Matrix Number) A)]{
To normalize a matrix is to scale it, such that the result has norm 1.
The call @racket[(matrix-normalize M p fail)] normalizes @racket[M] with respect to
the @racket[Lp]-norm. If the matrix @racket[M] has norm 0, the result of calling
the thunk @racket[fail] is returned.
If no fail-thunk is given, an argument-error exception is raised.
If no @racket[p] the L2-norm (Euclidean norm) is used.
@examples[#:eval untyped-eval
(require racket/math) ; for radians->degrees
@defproc[(matrix-normalize [M (Matrix Number)] [p Real 2] [fail (-> F) (λ () (error ...))])
(U F (Matrix Number))]{
Normalizes @racket[M] with respect to the L@subscript{@racket[p]} norm.
@examples[#:eval typed-eval
(matrix-normalize (col-matrix [1 1]))
(matrix-normalize (col-matrix [1 1]) 1)
(matrix-normalize (col-matrix [1 1]) +inf.0)]
The result of applying the @tech{failure thunk} @racket[fail] is returned if @racket[M]'s norm is
zero.
}
@deftogether[(@defproc[(matrix-normalize-rows [M (matrix Number)] [p Real 2] [fail (-> A) raise-argument-error]) (Matrix Number)]
@defproc[(matrix-normalize-cols [M (matrix Number)] [p Real 2] [fail (-> A) raise-argument-error]) (Matrix Number)])]{
@deftogether[(@defproc[(matrix-normalize-rows [M (Matrix Number)]
[p Real 2]
[fail (-> F) (λ () (error ...))]) (Matrix Number)]
@defproc[(matrix-normalize-cols [M (Matrix Number)]
[p Real 2]
[fail (-> F) (λ () (error ...))]) (Matrix Number)])]{
As @racket[matrix-normalize] but each row or column is normalized separately.
The result is this a matrix with unit vectors as rows or columns.
@examples[#:eval untyped-eval
(require racket/math) ; for radians->degrees
The result is a matrix with unit vectors as rows or columns.
@examples[#:eval typed-eval
(matrix-normalize-rows (matrix [[1 2] [2 4]]))
(matrix-normalize-cols (matrix [[1 2] [2 4]]))]
The result of applying the @tech{failure thunk} @racket[fail] is returned if the norm of any row
or column in @racket[M] is zero.
}
@deftogether[(@defproc[(matrix-rows-orthogonal? [M (Matrix Number)] [eps Real (* 10 epsilon.0)]) Boolean]
@defproc[(matrix-cols-orthogonal? [M (Matrix Number)] [eps Real (* 10 epsilon.0)]) Boolean])]{
@deftogether[(@defproc[(matrix-rows-orthogonal? [M (Matrix Number)]
[eps Real (* 10 epsilon.0)])
Boolean]
@defproc[(matrix-cols-orthogonal? [M (Matrix Number)]
[eps Real (* 10 epsilon.0)])
Boolean])]{
Returns @racket[#t] if the rows or columns of @racket[M]
are very close of being orthogonal (by default a few epsilons).
@examples[#:eval untyped-eval
(require racket/math) ; for radians->degrees
(matrix-rows-orthogonal? (matrix [[1 1] [-1 1]]))
(matrix-cols-orthogonal? (matrix [[1 1] [-1 1]]))]
}
@;{==================================================================================================}
@section[#:tag "matrix:solve"]{Solving Systems of Equations}
@defthing[matrix-solve Procedure]{}
@defproc[(matrix-solve [M (Matrix Number)] [B (Matrix Number)] [fail (-> F) (λ () (error ...))])
(U F (Matrix Number))]{
Returns the matrix @racket[X] for which @racket[(matrix* M X)] is @racket[B].
@racket[M] and @racket[B] must have the same number of rows.
@defthing[matrix-inverse Procedure]{}
It is typical for @racket[B] (and thus @racket[X]) to be a @tech{column matrix}, but not required.
If @racket[B] is not a column matrix, @racket[matrix-solve] solves for all the columns in @racket[B]
simultaneously.
@defthing[matrix-invertible? Procedure]{}
@examples[#:eval typed-eval
(define M (matrix [[7 5] [3 -2]]))
(define B0 (col-matrix [3 22]))
(define B1 (col-matrix [19 4]))
(matrix-solve M B0)
(matrix* M (col-matrix [4 -5]))
(matrix-solve M B1)
(matrix-cols (matrix-solve M (matrix-augment (list B0 B1))))]
@defthing[matrix-determinant Procedure]{}
@racket[matrix-solve] does not solve overconstrained or underconstrained systems, meaning that
@racket[M] must be invertible.
If @racket[M] is not invertible, the result of applying the @tech{failure thunk} @racket[fail] is
returned.
@racket[matrix-solve] is implemented using @racket[matrix-gauss-elim] to preserve exactness in its
output, with partial pivoting for greater numerical stability when @racket[M] is not exact.
See @racket[vandermonde-matrix] for an example that uses @racket[matrix-solve] to compute Legendre
polynomials.
}
@define[inverse-url]{http://en.wikipedia.org/wiki/Invertible_matrix}
@margin-note*{@hyperlink[inverse-url]{Wikipedia: Invertible Matrix}}
@defproc[(matrix-inverse [M (Matrix Number)] [fail (-> F) (λ () (error ...))]) (U F (Matrix Number))]{
Returns the inverse of @racket[M] if it exists; otherwise returns the result of applying the
@tech{failure thunk} @racket[fail].
@examples[#:eval typed-eval
(matrix-inverse (identity-matrix 3))
(matrix-inverse (matrix [[7 5] [3 -2]]))
(matrix-inverse (matrix [[1 2] [10 20]]))
(matrix-inverse (matrix [[1 2] [10 20]]) (λ () #f))]
}
@defproc[(matrix-invertible? [M (Matrix Number)]) Boolean]{
Returns @racket[#t] when @racket[M] is a @racket[square-matrix?] and @racket[(matrix-determinant M)]
is nonzero.
}
@margin-note*{@hyperlink["http://en.wikipedia.org/wiki/Determinant"]{Wikipedia: Determinant}}
@defproc[(matrix-determinant [M (Matrix Number)]) Number]{
Returns the determinant of @racket[M], which must be a @racket[square-matrix?].
@examples[#:eval typed-eval
(matrix-determinant (diagonal-matrix '(1 2 3 4)))
(* 1 2 3 4)
(matrix-determinant (matrix [[1 2] [10 20]]))
(matrix-determinant (col-matrix [1 2]))]
}
@;{==================================================================================================}
@ -609,11 +759,87 @@ are very close of being orthogonal (by default a few epsilons).
@section[#:tag "matrix:row-alg"]{Row-Based Algorithms}
@defthing[matrix-gauss-elim Procedure]{}
@define[gauss-url]{http://en.wikipedia.org/wiki/Gaussian_elimination}
@define[gauss-jordan-url]{http://en.wikipedia.org/wiki/Gauss%E2%80%93Jordan_elimination}
@defthing[matrix-row-echelon Procedure]{}
@margin-note*{@hyperlink[gauss-url]{Wikipedia: Gaussian elimination}}
@defproc[(matrix-gauss-elim [M (Matrix Number)]
[jordan? Any #f]
[unitize-pivot? Any #f]
[pivoting (U 'first 'partial) 'partial])
(Values (Matrix Number) (Listof Index))]{
@margin-note*{@hyperlink[gauss-jordan-url]{Wikipedia: Gauss-Jordan elimination}}
Implements Gaussian elimination or Gauss-Jordan elimination.
@defthing[matrix-lu Procedure]{}
If @racket[jordan?] is true, row operations are done both above and below the pivot.
If @racket[unitize-pivot?] is true, the pivot's row is scaled so that the pivot value is @racket[1].
When both are true, the algorithm is called Gauss-Jordan elimination, and the result matrix is in
@italic{reduced} row echelon form.
If @racket[pivoting] is @racket['first], the first nonzero entry in the current column is used as the
pivot. If @racket[pivoting] is @racket['partial], the largest-magnitude nonzero entry is used, which
improves numerical stability on average when @racket[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 @racket[matrix-row-echelon] for examples.
}
@define[row-echelon-url]{http://en.wikipedia.org/wiki/Row_echelon_form}
@margin-note*{@hyperlink[row-echelon-url]{Wikipedia: Row echelon form}}
@defproc[(matrix-row-echelon [M (Matrix Number)]
[jordan? Any #f]
[unitize-pivot? Any #f]
[pivoting (U 'first 'partial) 'partial])
(Matrix Number)]{
Like @racket[matrix-gauss-elim], but returns only the result of Gaussian elimination.
@examples[#:eval typed-eval
(define M (matrix [[2 1 -1] [-3 -1 2] [-2 1 2]]))
(matrix-row-echelon M)
(matrix-row-echelon M #t)
(matrix-identity? (matrix-row-echelon M #t #t))]
The last example shows that @racket[M] is invertible.
Using @racket[matrix-row-echelon] to solve a system of linear equations (without checking for
failure):
@interaction[#:eval typed-eval
(define B (col-matrix [8 -11 -3]))
(define MB (matrix-augment (list M B)))
(matrix-col (matrix-row-echelon MB #t #t) 3)
(matrix-solve M B)]
Using @racket[matrix-row-echelon] to invert a matrix (also without checking for failure):
@interaction[#:eval typed-eval
(define MI (matrix-augment (list M (identity-matrix 3))))
(submatrix (matrix-row-echelon MI #t #t) (::) (:: 3 #f))
(matrix-inverse M)]
}
@define[lu-url]{http://en.wikipedia.org/wiki/LU_decomposition}
@margin-note*{@hyperlink[lu-url]{Wikipedia: LU decomposition}}
@defproc[(matrix-lu [M (Matrix Number)] [fail (-> F) (λ () (error ...))])
(Values (U F (Matrix Number)) (Matrix Number))]{
Returns the LU decomposition of @racket[M] (which must be a @racket[square-matrix?]) if one exists.
An LU decomposition exists if @racket[M] can be put in row-echelon form without swapping rows.
Because @racket[matrix-lu] returns a @italic{unit} lower-triangular matrix (i.e. a lower-triangular
matrix with only ones on the diagonal), the decomposition is unique if it exists.
@examples[#:eval typed-eval
(define-values (L U)
(matrix-lu (matrix [[4 3] [6 3]])))
(values L U)
(matrix* L U)]
If @racket[M] does not have an LU decomposition, the first result is the result of applying the
@tech{failure thunk} @racket[fail], and the second result is the original argument @racket[M]:
@interaction[#:eval typed-eval
(matrix-lu (matrix [[0 1] [1 1]]))
(matrix-lu (matrix [[0 1] [1 1]]) (λ () #f))]
}
@;{==================================================================================================}
@ -621,57 +847,143 @@ are very close of being orthogonal (by default a few epsilons).
@section[#:tag "matrix:ortho-alg"]{Orthogonal Algorithms}
@defthing[matrix-gram-schmidt Procedure]{}
@define[gram-schmidt-url]{http://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process}
@define[reortho-pdf]{http://www.cerfacs.fr/algor/reports/2002/TR_PA_02_33.pdf}
@defthing[matrix-basis-extension Procedure]{}
@margin-note*{@hyperlink[gram-schmidt-url]{Wikipedia: Gram-Schmidt process}}
@defproc[(matrix-gram-schmidt [M (Matrix Number)] [normalize? Any #f] [start-col Integer 0])
(Array Number)]{
Returns an array whose columns are orthogonal and span the same subspace as @racket[M]'s columns.
The number of columns in the result is the rank of @racket[M].
If @racket[normalize?] is true, the columns are also normalized.
@margin-note{@hyperlink["http://en.wikipedia.org/wiki/QR_decomposition"]{Wikipedia: QR decomposition}}
@defproc*[([(matrix-qr [M (Matrix Real)]) (Values (Matrix Real) (Matrix Real))]
[(matrix-qr [M (Matrix Real)] [full Any]) (Values (Matrix Real) (Matrix Real))]
[(matrix-qr [M (Matrix Number)]) (Values (Matrix Number) (Matrix Number))]
[(matrix-qr [M (Matrix Number)] [full Any]) (Values (Matrix Number) (Matrix Number))])]{
@examples[#:eval typed-eval
(define M
(matrix [[12 -51 4]
[ 6 167 -68]
[-4 24 -41]]))
(matrix-gram-schmidt M)
(matrix-gram-schmidt M #t)
(matrix-cols-orthogonal? (matrix-gram-schmidt M))
(matrix-orthonormal? (matrix-gram-schmidt M #t))]
Examples with rank-deficient matrices:
@interaction[#:eval typed-eval
(matrix-gram-schmidt (matrix [[ 1 -2 1]
[-2 4 9]
[ 3 -6 7]]))
(matrix-gram-schmidt (make-matrix 3 3 0))]
When @racket[start-col] is positive, the Gram-Schmidt process is begun on column @racket[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
@racket[matrix-basis-extension].
@margin-note*{* @bold{On the round-off error analysis of the Gram-Schmidt algorithm with
reorthogonalization.},
Luc Giraud, Julien Langou and Miroslav Rozloznik. 2002. @hyperlink[reortho-pdf]{(PDF)}}
While Gram-Schmidt with inexact matrices is known to be unstable, using it twice tends to remove
instabilities:*
@interaction[#:eval typed-eval
(define M (matrix [[0.70000 0.70711]
[0.70001 0.70711]]))
(matrix-orthonormal?
(matrix-gram-schmidt M #t))
(matrix-orthonormal?
(matrix-gram-schmidt (matrix-gram-schmidt M) #t))]
}
@defproc[(matrix-basis-extension [M (Matrix Number)]) (Array Number)]{
Returns additional orthogonal columns which, if augmented with @racket[M], would result in an
orthogonal matrix of full rank. If @racket[M]'s columns are normalized, the result's columns are
normalized.
}
@define[qr-url]{http://en.wikipedia.org/wiki/QR_decomposition}
@margin-note*{@hyperlink["http://en.wikipedia.org/wiki/QR_decomposition"]{Wikipedia: QR decomposition}}
@defproc*[([(matrix-qr [M (Matrix Number)]) (Values (Matrix Number) (Matrix Number))]
[(matrix-qr [M (Matrix Number)] [full? Any]) (Values (Matrix Number) (Matrix Number))])]{
Computes a QR-decomposition of the matrix @racket[M]. The values returned are
the matrices @racket[Q] and @racket[R]. If @racket[full] is false, then
the matrices @racket[Q] and @racket[R]. If @racket[full?] is @racket[#f], then
a reduced decomposition is returned, otherwise a full decomposition is returned.
@margin-note{An @italic{orthonormal} matrix has columns which are orthooginal, unit vectors.}
@margin-note*{An @italic{orthonormal} matrix has columns which are orthoginal, unit vectors.}
The (full) decomposition of a square matrix consists of two matrices:
a orthogonal matrix @racket[Q] and an upper triangular matrix @racket[R],
such that @racket[QR = M].
For tall non-square matrices @racket[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 @hyperlink["http://en.wikipedia.org/wiki/QR_decomposition"]{QR decomposition}
See the Wikipedia entry on @hyperlink[qr-url]{QR decomposition}
for more details.
@examples[#:eval typed-eval
(define M
(matrix [[12 -51 4]
[ 6 167 -68]
[-4 24 -41]]))
(define-values (Q R) (matrix-qr M))
(values Q R)
(matrix= M (matrix* Q R))]
The difference between full and reduced decompositions:
@interaction[#:eval typed-eval
(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 Q2 R2)
(matrix= M (matrix* Q1 R1))
(matrix= M (matrix* Q2 R2))]
The decomposition @racket[M = QR] is useful for solving the equation @racket[Mx=v].
Since the inverse of Q is simply the transpose of Q,
@racket[Mx=v <=> QRx=v <=> Rx = Q^T v].
And since @racket[R] is upper triangular, the system can be solved by back substitution.
The algorithm used is Gram-Schmidt with reorthogonalization.
See the paper @hyperlink["http://www.cerfacs.fr/algor/reports/2002/TR_PA_02_33.pdf"]{On the round-off error analysis of the Gram-Schmidt algorithm with reorthogonalization.}
by Luc Giraud, Julien Langou, and, Miroslav Rozloznik.
}
}
@;{==================================================================================================}
@define[op-norm-url]{http://en.wikipedia.org/wiki/Matrix_norm#Induced_norm}
@section[#:tag "matrix:op-norm"]{Operator Norms and Comparing Matrices}
@secref{matrix:inner} describes functions that deal with matrices as vectors in
an inner product space. This section describes functions that deal with matrices as
@italic{linear operators}, or as functions from column matrices to column matrices.
@margin-note*{@hyperlink[op-norm-url]{Wikipedia: Induced norm}}
In this setting, a norm is the largest relative change in magnitude an operator (i.e. matrix) can
effect on a column matrix, where ``magnitude'' is defined by a vector norm.
(See the Wikipedia article linked to in the margin for a formal definition.)
Matrix norms that are defined in terms of a vector norm are called @italic{induced norms}, or
@deftech{operator norms}.
@defproc[(matrix-op-1norm [M (Matrix Number)]) Nonnegative-Real]{
TODO: describe
The @tech{operator norm} induced by the vector norm @racket[matrix-1norm].
When M is a column matrix, @racket[(matrix-op-1norm M)] is equivalent to @racket[(matrix-1norm M)].
}
@defproc[(matrix-op-2norm [M (Matrix Number)]) Nonnegative-Real]{
TODO: describe (spectral norm)
The @tech{operator norm} induced by the vector norm @racket[matrix-2norm].
This function is currently undefined because a required algorithm (singular value decomposition or
eigendecomposition) is not yet implemented in @racketmodname[math/matrix].
When M is a column matrix, @racket[(matrix-op-2norm M)] is equivalent to @racket[(matrix-2norm M)].
}
@defproc[(matrix-op-inf-norm [M (Matrix Number)]) Nonnegative-Real]{
TODO: describe
The @tech{operator norm} induced by the vector norm @racket[matrix-inf-norm].
When M is a column matrix, @racket[(matrix-op-inf-norm M)] is equivalent to
@racket[(matrix-inf-norm M)].
@ -680,16 +992,22 @@ When M is a column matrix, @racket[(matrix-op-inf-norm M)] is equivalent to
@defproc[(matrix-basis-cos-angle [M0 (Matrix Number)] [M1 (Matrix Number)]) Number]{
Returns the cosine of the angle between the two subspaces spanned by @racket[M0] and @racket[M1].
This function is currently undefined because a required algorithm (singular value decomposition or
eigendecomposition) is not yet implemented in @racketmodname[math/matrix].
When @racket[M0] and @racket[M1] are column matrices, @racket[(matrix-basis-cos-angle M0 M1)] is
equivalent to @racket[(matrix-cos-angle M0 M1)].
}
@defproc[(matrix-basis-angle [M0 (Matrix Number)] [M1 (Matrix Number)]) Number]{
Equivalent to @racket[(acos (matrix-basis-cos-angle M0 M1))].
The function is currently undefined because @racket[matrix-basis-cos-angle] is currently undefined.
}
@defparam[matrix-error-norm norm ((Matrix Number) -> Nonnegative-Real)]{
The norm used by @racket[matrix-relative-error] and @racket[matrix-absolute-error].
The default value is @racket[matrix-op-inf-norm].
Besides being a true norm, @racket[norm] should also be @deftech{submultiplicative}:
@racketblock[(norm (matrix* M0 M1)) <= (* (norm A) (norm B))]
@ -740,7 +1058,7 @@ Equivalent to
@defproc[(matrix-orthonormal? [M (Matrix Number)] [eps Real (* 10 epsilon.0)]) Boolean]{
Returns @racket[#t] when @racket[M] is very close to being orthonormal; that is, when
@racket[(matrix* M (matrix-hermitian M))] is very close to an identity matrix.
In fact, @racket[(matrix-orthonormal? M eps)] is equivalent to
Equivalent to
@racketblock[(matrix-identity? (matrix* M (matrix-hermitian M)) eps)]
}

View File

@ -17,17 +17,18 @@ for working with numbers and collections of numbers. These include
@item{@racketmodname[math/bigfloat]: Arbitrary-precision floating-point functions}
@item{@racketmodname[math/number-theory]: Number-theoretic functions}
@item{@racketmodname[math/array]: Functional arrays for operating on large rectangular data sets}
@item{@racketmodname[math/matrix]: Linear algebra functions for arrays @bold{(currently being documented)}}
@item{@racketmodname[math/matrix]: Linear algebra functions for arrays}
@item{@racketmodname[math/distributions]: Probability distributions}
@item{@racketmodname[math/statistics]: Statistical functions}
]
With this library, we hope to support a wide variety of applied mathematics in
Racket, including simulation, statistical inference, signal processing, and
combinatorics. If you find it lacking for your variety of mathematics,
please contact us or post to one of the
@hyperlink["http://racket-lang.org/community.html"]{mailing lists}
to make suggestions or submit patches.
combinatorics. If you find it lacking for your variety of mathematics, please
@itemlist[@item{Visit the @hyperlink["https://github.com/plt/racket/wiki/Math-Library-Features"]{Math Library Features wiki page}
to see what is planned.}
@item{Contact us or post to one of the @hyperlink["http://racket-lang.org/community.html"]{mailing lists} to
make suggestions or submit patches.}]
@bold{This is a Typed Racket library.} It is most efficient to use it in Typed
Racket, so that contracts are checked statically. However, almost all of it can

View File

@ -96,7 +96,7 @@
(check-true (array-strict? (matrix-col M i))))
(for*: ([M (list nonstrict-2x2-arr strict-2x2-arr)]
[spec (list '(0) 0)])
[spec (list '(0) (:: #f #f 2))])
(check-true (array-strict? (submatrix M (::) spec))))
)
@ -153,7 +153,7 @@
(check-false (array-strict? R))))
(for*: ([M (list nonstrict-2x2-arr strict-2x2-arr)]
[spec (list '(0) 0)])
[spec (list '(0) (:: #f #f 2))])
(check-false (array-strict? (submatrix M (::) spec))))
(for*: ([M (list nonstrict-2x2-arr strict-2x2-arr)]