Added QR decomposition tests (mostly a pretty thorough randomized test);

fixed error in `matrix-cols-orthogonal?' revealed by testing

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

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))
+(: pairwise-orthogonal? ((Listof (Matrix Number)) Nonnegative-Real -> Boolean))
-(define (rows-orthogonal? M eps)
+(define (pairwise-orthogonal? Ms eps)
-  (define rows (matrix->vector* M))
+  (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)]))

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"