Acos in matrix angle

math-lib/math/private/matrix/matrix-basic.rkt

@@ -309,15 +309,73 @@
                             ((Matrix Real) (Matrix Real) -> Real)
                             ((Matrix Float-Complex) (Matrix Float-Complex) -> Float-Complex)
                             ((Matrix Number) (Matrix Number) -> Number)))
-(define (matrix-cos-angle M N)
-  (/ (matrix-dot M N) (* (matrix-2norm M) (matrix-2norm N))))
+(define matrix-cos-angle
+  (let ()
+    (define (mag² [x : Number]) (if (real? x) (sqr x) (+ (sqr (real-part x)) (sqr (imag-part x)))))
+    (define (inf=>1 [x : Real]) : Flonum (if (eq? x +inf.0) 1. (if (eq? x -inf.0) -1. 0.)))
+    (: cinf=>1 (case-> (-> Real Flonum)
+                       (-> Float-Complex Float-Complex)
+                       (-> Number Number)))
+    (define (cinf=>1 x)
+      (if (real? x)
+          (inf=>1 x)
+          (make-flrectangular (inf=>1 (real-part x))
+                              (inf=>1 (imag-part x)))))
+    (: unit-bound (case-> (-> Flonum Flonum)
+                          (-> Real Real)))
+    (define (unit-bound x)
+      (if (flonum? x)
+          (flmin (flmax -1. x) 1.)
+          (min (max -1 x) 1)))
+    (: inner (case-> ((Matrix Flonum) (Matrix Flonum) -> Flonum)
+                     ((Matrix Real) (Matrix Real) -> Real)
+                     ((Matrix Float-Complex) (Matrix Float-Complex) -> Float-Complex)
+                     ((Matrix Number) (Matrix Number) -> Number)))
+    (define (inner A* B*)
+      (define nA (sqrt (array-all-sum (inline-array-map mag² A*))))
+      (define nB (sqrt (array-all-sum (inline-array-map mag² B*))))
+
+      (define result (/ (matrix-dot A* B*) (* nA nB)))
+
+      (if (real? result)
+          (unit-bound result)
+          (make-rectangular (unit-bound (real-part result))
+                            (unit-bound (imag-part result)))))
+    (λ (A B)
+      (define mA (array-strict (inline-array-map nonstupid-magnitude A)))
+      (define mB (array-strict (inline-array-map nonstupid-magnitude B)))
+      (define mxA (array-all-max mA))
+      (define mxB (array-all-max mB))
+
+      (cond
+        [(and (rational? mxA) (positive? mxA)
+              (rational? mxB) (positive? mxB))
+         (define A* (inline-array-map (λ (x) (/ x mxA)) A))
+         (define B* (inline-array-map (λ (x) (/ x mxB)) B))
+
+         (inner A* B*)]
+        [(or (nan? mxA) (nan? mxB) (= 0 mxA) (= 0 mxB))
+         (/ (matrix-dot A B) (* mxA mxB))]
+        [else
+         (define A* (if (rational? mxA)
+                        (inline-array-map (λ (x) (/ x mxA)) A)
+                        (array-strict (inline-array-map cinf=>1 A))))
+         (define B* (if (rational? mxB)
+                        (inline-array-map (λ (x) (/ x mxB)) B)
+                        (array-strict (inline-array-map cinf=>1 B))))
+         (inner A* B*)]))))
 
 (: matrix-angle (case-> ((Matrix Flonum) (Matrix Flonum) -> Flonum)
                         ((Matrix Real) (Matrix Real) -> Real)
                         ((Matrix Float-Complex) (Matrix Float-Complex) -> Float-Complex)
                         ((Matrix Number) (Matrix Number) -> Number)))
-(define (matrix-angle M N)
-  (acos (matrix-cos-angle M N)))
+(define (matrix-angle A B)
+  (define ca (matrix-cos-angle A B))
+  (cond
+    [(flonum? ca) (flacos ca)]
+    [(real? ca)
+     (if (eq? ca 1) 0 (flacos (fl ca)))]
+    [else (acos ca)]))
 
 (: matrix-normalize
    (All (A) (case-> ((Matrix Flonum)             -> (Matrix Flonum))

math-lib/math/private/matrix/matrix-operator-norm.rkt

@@ -79,9 +79,7 @@ See "How to Measure Errors" in the LAPACK manual for more details:
   ;(matrix-min-singular-value (matrix* (matrix-hermitian M) R))
   (error 'unimplemented))
 
-(: matrix-basis-angle (case-> ((Matrix Flonum) (Matrix Flonum) -> Flonum)
-                              ((Matrix Real) (Matrix Real) -> Real)
-                              ((Matrix Float-Complex) (Matrix Float-Complex) -> Float-Complex)
+(: matrix-basis-angle (case-> ((Matrix Float-Complex) (Matrix Float-Complex) -> Float-Complex)
                               ((Matrix Number) (Matrix Number) -> Number)))
 ;; Returns the angle between the two subspaces spanned by the two given sets of column vectors
 (define (matrix-basis-angle M R)

math-test/math/tests/matrix-tests.rkt

@@ -661,7 +661,39 @@
   (check-exn exn:fail:contract? (λ () (matrix-dot a (matrix [[1]]))))
   (check-exn exn:fail:contract? (λ () (matrix-dot (matrix [[1]]) a))))
 
-;; TODO: matrix-angle
+;; matrix-angle
+;; guarantee that matrix-cos-angle on real is in [-1 .. 1] except for zero?
+(for: ([i (in-range 10)])
+      (define n (+ 1 (random 3)))
+      (define m (+ 1 (random 4)))
+      (define A (random-matrix n m))
+      (define B (random-matrix n m))
+      (define ca (matrix-cos-angle A B))
+      (check-true (or (and (rational? ca) (<= ca 1.))
+                      (matrix-zero? A) (matrix-zero? B))))
+(let ([A (matrix [[-1.65e+145 -9.55e+152]])]
+      [B (matrix [[-2.05e-130 -1.88e-122]])])
+  (check-true (<= (abs (matrix-cos-angle A B)) 1.)))
+;; check matrix-cos-angle doesn't overflow easily
+(let ([A (matrix [[-1.20e-21+8.16e+173i -3.18e+280-4.15e+275i]])]
+      [B (matrix [[-2.51e-270-4.38e-68i -1.65e+232+3.06e+227i]])])
+  (define ca (matrix-cos-angle A B))
+  ;; 0.9999999995008538+3.159576900273938e-05i
+  (check-true (<= (magnitude ca) (flnext 1.))))
+(let ([A (matrix-scale (matrix [[1 1]]) 1e270+3.i)]
+      [B (matrix-scale (matrix [[9 0]]) 1e270+3.i)])
+  (define ca (matrix-cos-angle A B))
+  ;; 0.9999999995008538+3.159576900273938e-05i
+  (check-true (<= (magnitude (- ca (make-rectangular (sqrt 1/2) 0.))) epsilon.0)))
+;; check matrix-cos-angle treats axes with inf as principal axes
+(check-equal? (matrix-cos-angle (matrix [[1 5] [2 9]])
+                                (matrix [[2 8] [3 -inf.0]]))
+              (matrix-cos-angle (matrix [[1 5] [2 9]])
+                                (matrix [[0. 0.] [0. -1.0]])))
+(check-equal? (matrix-cos-angle (matrix [[2 9]])
+                                (matrix [[3 -8+inf.0i]]))
+              (matrix-cos-angle (matrix [[2 9]])
+                                (matrix [[0. -0+1.0i]])))
 
 ;; TODO: matrix-normalize