diff --git a/collects/math/matrix.rkt b/collects/math/matrix.rkt index e3829832eb..b238016fca 100644 --- a/collects/math/matrix.rkt +++ b/collects/math/matrix.rkt @@ -18,6 +18,7 @@ vandermonde-matrix) (except-in "private/matrix/matrix-basic.rkt" matrix-dot + matrix-cos-angle matrix-angle matrix-normalize matrix-conjugate @@ -70,6 +71,8 @@ [matrix-dot (case-> ((Matrix Number) -> Nonnegative-Real) ((Matrix Number) (Matrix Number) -> Number))] + [matrix-cos-angle + ((Matrix Number) (Matrix Number) -> Number)] [matrix-angle ((Matrix Number) (Matrix Number) -> Number)] [matrix-normalize diff --git a/collects/math/private/matrix/matrix-basic.rkt b/collects/math/private/matrix/matrix-basic.rkt index 3fee71a986..8611455276 100644 --- a/collects/math/private/matrix/matrix-basic.rkt +++ b/collects/math/private/matrix/matrix-basic.rkt @@ -40,6 +40,7 @@ matrix-inf-norm matrix-norm matrix-dot + matrix-cos-angle matrix-angle matrix-normalize ;; Simple operators @@ -240,10 +241,15 @@ (λ: ([js : Indexes]) (* (aproc js) (conjugate (bproc js))))))])) +(: matrix-cos-angle (case-> ((Matrix Real) (Matrix Real) -> Real) + ((Matrix Number) (Matrix Number) -> Number))) +(define (matrix-cos-angle M N) + (/ (matrix-dot M N) (* (matrix-2norm M) (matrix-2norm N)))) + (: matrix-angle (case-> ((Matrix Real) (Matrix Real) -> Real) ((Matrix Number) (Matrix Number) -> Number))) (define (matrix-angle M N) - (acos (/ (matrix-dot M N) (* (matrix-2norm M) (matrix-2norm N))))) + (acos (matrix-cos-angle M N))) (: matrix-normalize (All (A) (case-> ((Matrix Real) -> (Matrix Real)) @@ -361,7 +367,7 @@ (matrix-map-cols (make-matrix-normalize p) M fail)])) ;; =================================================================================================== -;; Robust predicates using entrywise norms +;; Approximate predicates using entrywise norms (: matrix-zero? (case-> ((Matrix Number) -> Boolean) ((Matrix Number) Real -> Boolean))) @@ -369,26 +375,26 @@ (cond [(negative? eps) (raise-argument-error 'matrix-zero? "Nonnegative-Real" 1 M eps)] [else (<= (matrix-norm M +inf.0) eps)])) -(: rows-orthogonal? ((Matrix Number) Nonnegative-Real -> Boolean)) -(define (rows-orthogonal? M eps) - (define rows (matrix->vector* M)) +(: pairwise-orthogonal? ((Listof (Matrix Number)) Nonnegative-Real -> Boolean)) +(define (pairwise-orthogonal? Ms eps) + (define rows (list->vector Ms)) (define m (vector-length rows)) (let/ec: return : Boolean (for*: ([i0 (in-range m)] [i1 (in-range (fx+ i0 1) m)]) (define r0 (unsafe-vector-ref rows i0)) (define r1 (unsafe-vector-ref rows i1)) - (when ((sqrt (magnitude (vector-dot r0 r1))) . >= . eps) (return #f))) + (when ((magnitude (matrix-cos-angle r0 r1)) . >= . eps) (return #f))) #t)) (: matrix-rows-orthogonal? (case-> ((Matrix Number) -> Boolean) ((Matrix Number) Real -> Boolean))) (define (matrix-rows-orthogonal? M [eps (* 10 epsilon.0)]) (cond [(negative? eps) (raise-argument-error 'matrix-rows-orthogonal? "Nonnegative-Real" 1 M eps)] - [else (rows-orthogonal? M eps)])) + [else (pairwise-orthogonal? (matrix-rows M) eps)])) (: matrix-cols-orthogonal? (case-> ((Matrix Number) -> Boolean) ((Matrix Number) Real -> Boolean))) (define (matrix-cols-orthogonal? M [eps (* 10 epsilon.0)]) (cond [(negative? eps) (raise-argument-error 'matrix-cols-orthogonal? "Nonnegative-Real" 1 M eps)] - [else (rows-orthogonal? (matrix-transpose M) eps)])) + [else (pairwise-orthogonal? (matrix-cols M) eps)])) diff --git a/collects/math/tests/matrix-tests.rkt b/collects/math/tests/matrix-tests.rkt index 55340d7434..81edb48b33 100644 --- a/collects/math/tests/matrix-tests.rkt +++ b/collects/math/tests/matrix-tests.rkt @@ -930,14 +930,6 @@ [ 3/7 158/175 6/175] [-2/7 6/35 -33/35 ]])) -(check-equal? (matrix-gram-schmidt (matrix [[12 -51 4] - [ 6 167 -68] - [-4 24 -41]]) - #t) - (matrix [[ 6/7 -69/175 -58/175] - [ 3/7 158/175 6/175] - [-2/7 6/35 -33/35 ]])) - (check-equal? (matrix-gram-schmidt (matrix [[12 -51] [ 6 167] [-4 24]]) @@ -955,36 +947,67 @@ ;; =================================================================================================== ;; QR decomposition -(check-true (matrix-orthonormal? (matrix-q (index-array #(100 1))))) +(check-true (matrix-orthonormal? (matrix-q (index-array #(50 1))))) + +(let-values ([(Q R) (matrix-qr (matrix [[12 -51 4] + [ 6 167 -68] + [-4 24 -41]]))]) + (check-equal? Q (matrix [[ 6/7 -69/175 -58/175] + [ 3/7 158/175 6/175] + [-2/7 6/35 -33/35 ]])) + (check-equal? R (matrix [[14 21 -14] + [ 0 175 -70] + [ 0 0 35]]))) + +;; A particularly tricky test case used to demonstrate loss of orthogonality +;; QR has to generate a better Q than Gram-Schmidt alone (which fails this test) +(check-true (matrix-orthonormal? + (matrix-q (matrix [[0.70000 0.70711] + [0.70001 0.70711]])))) + +;; Fuzz test the heck out of it: 100 matrices, random shape, random entries, sometimes rank-deficient +(for: ([i (in-range 100)]) + (define m (+ 1 (random 10))) + (define n (+ 1 (random 10))) + (define M (random-matrix m n -3 4)) + ;; Full QR, real matrix + (let-values ([(Q R) (matrix-qr M #t)]) + (check-true (matrix-orthonormal? Q) + (format "M = ~a Q = ~a" M Q)) + (check-true (<= (matrix-relative-error (matrix* Q R) M) + (* 10 epsilon.0)))) + ;; Reduced QR, real matrix + (let-values ([(Q R) (matrix-qr M #f)]) + (check-true (matrix-cols-orthogonal? Q) + (format "M = ~a Q = ~a" M Q)) + (check-true (<= (matrix-relative-error (matrix* Q R) M) + (* 10 epsilon.0)))) + (define N (random-matrix m n -3 4)) + (define M+N (array-make-rectangular M N)) + ;; Full QR, complex matrix + (let-values ([(Q R) (matrix-qr M+N #t)]) + (check-true (matrix-orthonormal? Q) + (format "M+N = ~a Q = ~a" M+N Q)) + (check-true (<= (matrix-relative-error (matrix* Q R) M+N) + (* 10 epsilon.0)))) + ;; Reduced QR, complex matrix + (let-values ([(Q R) (matrix-qr M+N #f)]) + (check-true (matrix-cols-orthogonal? Q) + (format "M+N = ~a Q = ~a" M+N Q)) + (check-true (<= (matrix-relative-error (matrix* Q R) M+N) + (* 10 epsilon.0))))) + +(for: ([M (in-list nonmatrices)]) + (check-exn exn:fail:contract? (λ () (matrix-q M)))) #| ;; =================================================================================================== ;; Tests not yet converted to rackunit -;; A particularly tricky one used to demonstrate loss of orthogonality -(matrix-qr (matrix [[0.70000 0.70711] - [0.70001 0.70711]])) - (begin (begin "matrix-operations.rkt" - #; - (list 'column-dimension - (= (column-dimension #(1 2 3)) 3) - (= (column-dimension (vector->matrix 1 2 #(1 2))) 1)) - (let ([matrix: vector->matrix]) - (list 'column-dot - (= (column-dot (col-matrix [1 2]) (col-matrix [1 2])) 5) - (= (column-dot (col-matrix [1 2]) (col-matrix [3 4])) 11) - (= (column-dot (col-matrix [3 4]) (col-matrix [3 4])) 25) - (= (column-dot (col-matrix [1 2 3]) (col-matrix [4 5 6])) - (+ (* 1 4) (* 2 5) (* 3 6))) - (= (column-dot (col-matrix [+3i +4i]) (col-matrix [+3i +4i])) - 25))) - (let ([matrix: vector->matrix]) - (list 'column-norm - (= (column-norm (col-matrix [2 4])) (sqrt 20)))) (list 'column-project (equal? (column-project #(1 2 3) #(4 5 6)) (col-matrix [128/77 160/77 192/77])) (equal? (column-project (col-matrix [1 2 3]) (col-matrix [2 4 3])) @@ -1002,45 +1025,9 @@ (matrix-scale (col-matrix [-1 1 -1 1]) 1/2) (matrix-scale (col-matrix [ 1 -1 -1 1]) 1/2))) (col-matrix [2 3 2 3]))) - (list 'gram-schmidt-orthogonal - (equal? (gram-schmidt-orthogonal (list #(3 1) #(2 2))) - (list (col-matrix [3 1]) (col-matrix [-2/5 6/5])))) - (list 'vector-normalize - (equal? (column-normalize #(3 4)) - (col-matrix [3/5 4/5]))) - (list 'gram-schmidt-orthonormal - (equal? (gram-schmidt-orthonormal (ann '(#(3 1) #(2 2)) (Listof (Column Number)))) - (list (column-normalize #(3 1)) - (column-normalize #(-2/5 6/5))))) - (list 'projection-on-subspace (equal? (projection-on-subspace #(1 2 3) '(#(2 4 3))) (matrix-scale (col-matrix [2 4 3]) 19/29))) - (list 'unit-vector - (equal? (unit-column 4 1) (col-matrix [0 1 0 0]))) - (list 'matrix-qr - (let-values ([(Q R) (matrix-qr (matrix [[1 1] [0 1] [1 1]]))]) - (equal? (list Q R) - (list (matrix [[0.7071067811865475 0] - [0 1] - [0.7071067811865475 0]]) - (matrix [[1.414213562373095 1.414213562373095] - [0 1]])))) - (let () - (define A (matrix [[1 2 3 4] [1 2 4 5] [1 2 5 6] [1 2 6 7]])) - (define-values (Q R) (matrix-qr A)) - (equal? (list Q R) - (list - (vector->matrix - 4 4 ((inst vector Number) - 1/2 -0.6708203932499369 0.5477225575051662 -0.0 - 1/2 -0.22360679774997896 -0.7302967433402214 0.4082482904638629 - 1/2 0.22360679774997896 -0.18257418583505536 -0.8164965809277259 - 1/2 0.6708203932499369 0.3651483716701107 0.408248290463863)) - (vector->matrix - 4 4 ((inst vector Number) - 2 4 9 11 0 0.0 2.23606797749979 2.23606797749979 - 0 0 0.0 4.440892098500626e-16 0 0 0 0.0)))))) #; (begin "matrix-2d.rkt"