suzanne.soy/racket
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Slightly more `math/matrix'

Browse Source 
* Moved to-do list in "matrix-operations.rkt" to the wiki

* Added more mutating vector ops

* Added "matrix-basis.rkt" (unfinished)
This commit is contained in:
Neil Toronto 2012-12-27 17:28:17 -07:00
parent e55a31480e
commit 7ac8e1bbce
4 changed files with 207 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

112
collects/math/private/matrix/matrix-basis.rkt Normal file
View File

@ -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]]))
|#

16
collects/math/private/matrix/matrix-operations.rkt
View File

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

87
collects/math/private/vector/vector-mutate.rkt
View File

@ -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])))

10
collects/math/tests/matrix-tests.rkt
View File

@ -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