From 7ac8e1bbcef78f17c1839e1ed582c535852d53d0 Mon Sep 17 00:00:00 2001 From: Neil Toronto Date: Thu, 27 Dec 2012 17:28:17 -0700 Subject: [PATCH] Slightly more `math/matrix' * Moved to-do list in "matrix-operations.rkt" to the wiki * Added more mutating vector ops * Added "matrix-basis.rkt" (unfinished) --- collects/math/private/matrix/matrix-basis.rkt | 112 ++++++++++++++++++ .../math/private/matrix/matrix-operations.rkt | 16 +-- .../math/private/vector/vector-mutate.rkt | 87 +++++++++++++- collects/math/tests/matrix-tests.rkt | 10 +- 4 files changed, 207 insertions(+), 18 deletions(-) create mode 100644 collects/math/private/matrix/matrix-basis.rkt diff --git a/collects/math/private/matrix/matrix-basis.rkt b/collects/math/private/matrix/matrix-basis.rkt new file mode 100644 index 0000000000..75bc4288ea --- /dev/null +++ b/collects/math/private/matrix/matrix-basis.rkt @@ -0,0 +1,112 @@ +#lang typed/racket + +(require racket/fixnum + math/array + math/matrix + "matrix-column.rkt" + "utils.rkt" + "../unsafe.rkt" + "../vector/vector-mutate.rkt" + ) + +(: col-matrix-project1 (case-> ((Matrix Real) (Matrix Real) Any -> (U #f (Matrix Real))) + ((Matrix Number) (Matrix Number) Any -> (U #f (Matrix Number))))) +(define (col-matrix-project1 v b unit?) + (cond [unit? (matrix-scale b (matrix-dot v b))] + [else (define b.b (matrix-dot b b)) + (cond [(and (zero? b.b) (exact? b.b)) #f] + [else (matrix-scale b (/ (matrix-dot v b) b.b))])])) + +(: col-matrix-project + (All (A) (case-> ((Matrix Real) (Matrix Real) -> (Matrix Real)) + ((Matrix Real) (Matrix Real) Any -> (U A (Matrix Real))) + ((Matrix Real) (Matrix Real) Any (-> A) -> (U A (Matrix Real))) + ((Matrix Number) (Matrix Number) -> (Matrix Number)) + ((Matrix Number) (Matrix Number) Any -> (U A (Matrix Number))) + ((Matrix Number) (Matrix Number) Any (-> A) -> (U A (Matrix Number)))))) +(define col-matrix-project + (case-lambda + [(v B) (col-matrix-project v B #f)] + [(v B unit?) + (col-matrix-project + v B unit? + (λ () (error 'col-matrix-project "expected basis with nonzero column vectors; given ~e" B)))] + [(v B unit? fail) + (unless (col-matrix? v) (raise-argument-error 'col-matrix-project "col-matrix?" v)) + (define bs (matrix-cols (ensure-matrix 'col-matrix-project B))) + (define p (col-matrix-project1 v (first bs) unit?)) + (cond [p (let loop ([bs (rest bs)] [p p]) + (cond [(empty? bs) p] + [else (define q (col-matrix-project1 v (first bs) unit?)) + (if q (loop (rest bs) (matrix+ p q)) (fail))]))] + [else (fail)])])) + +(: find-nonzero-vector (case-> ((Vectorof (Vectorof Real)) -> (U #f Index)) + ((Vectorof (Vectorof Number)) -> (U #f Index)))) +(define (find-nonzero-vector vss) + (define n (vector-length vss)) + (cond [(= n 0) #f] + [else (let loop ([#{i : Nonnegative-Fixnum} 0]) + (cond [(i . fx< . n) + (define vs (unsafe-vector-ref vss i)) + (if (vector-zero? vs) (loop (fx+ i 1)) i)] + [else #f]))])) + +(: subtract-projections! + (case-> ((Vectorof (Vectorof Real)) Index Index (Vectorof Real) Any -> Void) + ((Vectorof (Vectorof Number)) Index Index (Vectorof Number) Any -> Void))) +(define (subtract-projections! cols n i ci unit?) + (let j-loop ([#{j : Nonnegative-Fixnum} (fx+ i 1)]) + (when (j . fx< . n) + (vector-sub-proj! (unsafe-vector-ref cols j) ci unit?) + (j-loop (fx+ j 1))))) + +(: matrix-gram-schmidt (All (A) (case-> ((Matrix Real) -> (Array Real)) + ((Matrix Real) Any -> (Array Real)) + ((Matrix Number) -> (Array Number)) + ((Matrix Number) Any -> (Array Number))))) +(define (matrix-gram-schmidt M [unit? #f]) + (define rows (matrix->vector* M)) + (define n (vector-length rows)) + (define i (find-nonzero-vector rows)) + (cond [i (define rowi (unsafe-vector-ref rows i)) + (subtract-projections! rows n i rowi #f) + (when unit? (vector-normalize! rowi)) + (let loop ([#{i : Nonnegative-Fixnum} (fx+ i 1)] [bs (list rowi)]) + (cond [(i . fx< . n) + (define rowi (unsafe-vector-ref rows i)) + (cond [(vector-zero? rowi) (loop (fx+ i 1) bs)] + [else (subtract-projections! rows n i rowi #f) + (when unit? (vector-normalize! rowi)) + (loop (fx+ i 1) (cons rowi bs))])] + [else + (vector*->matrix (list->vector (reverse bs)))]))] + [else + (make-array (vector 0 (matrix-num-cols M)) 0)])) +#| +(define a (col-matrix [1 2 1])) +(define b (col-matrix [1 -2 2])) + +(define basis + (gram-schmidt-orthogonal + (matrix-cols + (array #[#[2 1 0] #[2 2 1] #[0 2 0]])))) + +(column-project a b) +(col-matrix-project a b) + +(projection-on-orthogonal-basis a basis) +(col-matrix-project a (matrix-augment basis)) +(projection-on-orthonormal-basis a basis) +(col-matrix-project a (matrix-augment basis) 'orthonormal) + +(matrix-gram-schmidt + (matrix [[0 1 2] + [0 2 3] + [0 1 5]])) + +(matrix-gram-schmidt + (matrix [[5 1 2] + [2 2 3] + [-3 1 5]])) +|# diff --git a/collects/math/private/matrix/matrix-operations.rkt b/collects/math/private/matrix/matrix-operations.rkt index 1582d3f497..a2a994e394 100644 --- a/collects/math/private/matrix/matrix-operations.rkt +++ b/collects/math/private/matrix/matrix-operations.rkt @@ -16,14 +16,6 @@ "utils.rkt" (for-syntax racket)) -; TODO: -; 1. compute null space from QR factorization -; (better numerical stability than from Gauss elimination) -; 2. S+N decomposition -; 3. Linear least squares problems (data fitting) -; 4. Pseudo inverse -; 5. Eigenvalues and eigenvectors - (provide ;; Gaussian elimination matrix-gauss-elim @@ -33,7 +25,7 @@ matrix-nullity matrix-determinant matrix-determinant/row-reduction ; for testing - ;; Spaces (TODO: null space, row space, left null space) + ;; Spaces matrix-column-space ;; Solving matrix-invertible? @@ -141,8 +133,6 @@ (: matrix-rank : (Matrix Number) -> Index) ;; Returns the dimension of the column space (equiv. row space) of M (define (matrix-rank M) - ; TODO: Use QR or SVD instead for inexact matrices - ; See answer: http://scicomp.stackexchange.com/questions/1861/understanding-how-numpy-does-svd (define n (matrix-num-cols M)) (define-values (_ cols-without-pivot) (matrix-gauss-elim M)) (assert (- n (length cols-without-pivot)) index?)) @@ -333,7 +323,7 @@ (matrix-sum (map (λ: ([b : (Column Number)]) (column-project v (->col-matrix b))) bs)))) - + ; (projection-on-orthonormal-basis v bs) ; Project the vector v on the orthonormal basis vectors in bs. ; The basis bs must be either the column vectors of a matrix @@ -394,7 +384,7 @@ (: extend-span-to-basis : (Listof (Matrix Number)) Integer -> (Listof (Matrix Number))) -; Extend the basis in vs to with rdimensional basis +; Extend the basis in vs to r-dimensional basis (define (extend-span-to-basis vs r) (define-values (m n) (matrix-shape (car vs))) (: loop : (Listof (Matrix Number)) (Listof (Matrix Number)) Integer -> (Listof (Matrix Number))) diff --git a/collects/math/private/vector/vector-mutate.rkt b/collects/math/private/vector/vector-mutate.rkt index 4ca10c35e2..9e4356c419 100644 --- a/collects/math/private/vector/vector-mutate.rkt +++ b/collects/math/private/vector/vector-mutate.rkt @@ -1,11 +1,24 @@ #lang typed/racket/base (require racket/fixnum + racket/math math/private/unsafe) (provide vector-swap! vector-scale! - vector-scaled-add!) + vector-scaled-add! + vector-mag^2 + vector-dot + vector-normalize! + vector-sub-proj! + vector-zero! + vector-zero?) + +(: mag^2 (Number -> Nonnegative-Real)) +(define (mag^2 x) + (define y (* x (conjugate x))) + (cond [(and (real? y) (y . >= . 0)) y] + [else (error 'impossible)])) (: vector-swap! (All (A) ((Vectorof A) Integer Integer -> Void))) (define (vector-swap! vs i0 i1) @@ -43,5 +56,73 @@ (: vector-scaled-add! (case-> ((Vectorof Real) (Vectorof Real) Real -> Void) ((Vectorof Number) (Vectorof Number) Number -> Void))) -(define (vector-scaled-add! v0 v1 s) - (vector-generic-scaled-add! v0 v1 s + *)) +(define (vector-scaled-add! vs0 vs1 s) + (vector-generic-scaled-add! vs0 vs1 s + *)) + +(: vector-mag^2 (case-> ((Vectorof Real) -> Nonnegative-Real) + ((Vectorof Number) -> Nonnegative-Real))) +(define (vector-mag^2 vs) + (define n (vector-length vs)) + (let loop ([#{i : Nonnegative-Fixnum} 0] [#{s : Nonnegative-Real} 0]) + (if (i . fx>= . n) s (loop (fx+ i 1) (+ s (mag^2 (unsafe-vector-ref vs i))))))) + +(: vector-dot (case-> ((Vectorof Real) (Vectorof Real) -> Real) + ((Vectorof Number) (Vectorof Number) -> Number))) +(define (vector-dot vs0 vs1) + (define n (min (vector-length vs0) (vector-length vs1))) + (cond [(= n 0) 0] + [else + (define v0 (unsafe-vector-ref vs0 0)) + (define v1 (unsafe-vector-ref vs1 0)) + (let loop ([#{i : Nonnegative-Fixnum} 1] [s (* v0 (conjugate v1))]) + (cond [(i . fx< . n) + (define v0 (unsafe-vector-ref vs0 i)) + (define v1 (unsafe-vector-ref vs1 i)) + (loop (fx+ i 1) (+ s (* v0 (conjugate v1))))] + [else s]))])) + +(: vector-normalize! (case-> ((Vectorof Real) -> Nonnegative-Real) + ((Vectorof Number) -> Nonnegative-Real))) +(define (vector-normalize! vs) + (define n (vector-length vs)) + (define s (sqrt (vector-mag^2 vs))) + (unless (and (zero? s) (exact? s)) + (let loop ([#{i : Nonnegative-Fixnum} 0]) + (when (i . fx< . n) + (unsafe-vector-set! vs i (/ (unsafe-vector-ref vs i) s)) + (loop (fx+ i 1))))) + s) + +(: vector-sub-proj! (case-> ((Vectorof Real) (Vectorof Real) Any -> Nonnegative-Real) + ((Vectorof Number) (Vectorof Number) Any -> Nonnegative-Real))) +(define (vector-sub-proj! vs0 vs1 unit?) + (define n (min (vector-length vs0) (vector-length vs1))) + (define t (if unit? 1 (vector-mag^2 vs1))) + (unless (and (zero? t) (exact? t)) + (define s (/ (vector-dot vs0 vs1) t)) + (let loop ([#{i : Nonnegative-Fixnum} 0]) + (when (i . fx< . n) + (define v0 (unsafe-vector-ref vs0 i)) + (define v1 (unsafe-vector-ref vs1 i)) + (unsafe-vector-set! vs0 i (- v0 (* v1 s))) + (loop (fx+ i 1))))) + t) + +(: vector-zero! (case-> ((Vectorof Real) -> Void) + ((Vectorof Number) -> Void))) +(define (vector-zero! vs) + (define n (vector-length vs)) + (let loop ([#{i : Nonnegative-Fixnum} 0]) + (when (i . fx< . n) + (unsafe-vector-set! vs i 0) + (loop (fx+ i 1))))) + +(: vector-zero? (case-> ((Vectorof Real) -> Boolean) + ((Vectorof Number) -> Boolean))) +(define (vector-zero? vs) + (define n (vector-length vs)) + (let loop ([#{i : Nonnegative-Fixnum} 0]) + (cond [(i . fx>= . n) #t] + [(zero? (unsafe-vector-ref vs i)) + (loop (fx+ i 1))] + [else #f]))) diff --git a/collects/math/tests/matrix-tests.rkt b/collects/math/tests/matrix-tests.rkt index c58fba12dc..368b69a12e 100644 --- a/collects/math/tests/matrix-tests.rkt +++ b/collects/math/tests/matrix-tests.rkt @@ -20,7 +20,7 @@ (make-array #(0 1) 0) (make-array #(0 0) 0) (make-array #(1 1 1) 0))) -#| + ;; =================================================================================================== ;; Literal syntax @@ -796,7 +796,7 @@ (for: ([a (in-list nonmatrices)]) (check-exn exn:fail:contract? (λ () (matrix-inverse a)))) -|# + ;; =================================================================================================== ;; LU decomposition @@ -834,6 +834,12 @@ ;; =================================================================================================== ;; Tests not yet converted to rackunit +(matrix-gram-schmidt + (matrix [[2 1 2] + [2 2 3] + [5 1 5]]) + #t) + (begin (begin