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

Fixed major performance issue with matrix arithmetic; please merge to 5.3.2

Browse Source 
The fix consists of three parts:

1. Rewriting `inline-matrix*'. The material change here is that the
   expansion now contains only direct applications of `+' and `*'.
   TR's optimizer replaces them with `unsafe-fx+' and `unsafe-fx*',
   which keeps intermediate flonum values from being boxed.

2. Making the types of all functions that operate on (Matrix Number)
   values more precise. Now TR can prove that matrix operations preserve
   inexactness. For example, matrix-conjugate : (Matrix Flonum) ->
   (Matrix Flonum) and three other cases for Real, Float-Complex, and
   Number.

3. Changing the return types of some functions that used to return
   things like (Matrix (U A 0)). Now that we worry about preserving
   inexactness, we can't have `matrix-upper-triangle' always return a
   matrix that contains exact zeros. It now accepts an optional `zero'
   argument of type A.
This commit is contained in:
Neil Toronto 2013-01-21 21:55:22 -07:00
parent a0f910c3dc
commit f42cc6f14a
17 changed files with 606 additions and 219 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
collects/math/matrix.rkt
View File

@ -17,6 +17,10 @@
(except-in "private/matrix/matrix-constructors.rkt" (except-in "private/matrix/matrix-constructors.rkt"
vandermonde-matrix) vandermonde-matrix)
(except-in "private/matrix/matrix-basic.rkt" (except-in "private/matrix/matrix-basic.rkt"
matrix-1norm
matrix-2norm
matrix-inf-norm
matrix-norm
matrix-dot matrix-dot
matrix-cos-angle matrix-cos-angle
matrix-angle matrix-angle
@ -29,6 +33,9 @@
(except-in "private/matrix/matrix-subspace.rkt" (except-in "private/matrix/matrix-subspace.rkt"
matrix-col-space) matrix-col-space)
(except-in "private/matrix/matrix-operator-norm.rkt" (except-in "private/matrix/matrix-operator-norm.rkt"
matrix-op-1norm
matrix-op-2norm
matrix-op-inf-norm
matrix-basis-cos-angle matrix-basis-cos-angle
matrix-basis-angle) matrix-basis-angle)
;;"private/matrix/matrix-qr.rkt" ; all use require/untyped-contract ;;"private/matrix/matrix-qr.rkt" ; all use require/untyped-contract
@ -77,6 +84,11 @@
(require/untyped-contract (require/untyped-contract
(begin (require "private/matrix/matrix-types.rkt")) (begin (require "private/matrix/matrix-types.rkt"))
"private/matrix/matrix-basic.rkt" "private/matrix/matrix-basic.rkt"
[matrix-1norm ((Matrix Number) -> Nonnegative-Real)]
[matrix-2norm ((Matrix Number) -> Nonnegative-Real)]
[matrix-inf-norm ((Matrix Number) -> Nonnegative-Real)]
[matrix-norm (case-> ((Matrix Number) -> Nonnegative-Real)
((Matrix Number) Real -> Nonnegative-Real))]
[matrix-dot [matrix-dot
(case-> ((Matrix Number) -> Nonnegative-Real) (case-> ((Matrix Number) -> Nonnegative-Real)
((Matrix Number) (Matrix Number) -> Number))] ((Matrix Number) (Matrix Number) -> Number))]
@ -113,6 +125,9 @@
(require/untyped-contract (require/untyped-contract
(begin (require "private/matrix/matrix-types.rkt")) (begin (require "private/matrix/matrix-types.rkt"))
"private/matrix/matrix-operator-norm.rkt" "private/matrix/matrix-operator-norm.rkt"
[matrix-op-1norm ((Matrix Number) -> Nonnegative-Real)]
[matrix-op-2norm ((Matrix Number) -> Nonnegative-Real)]
[matrix-op-inf-norm ((Matrix Number) -> Nonnegative-Real)]
[matrix-basis-cos-angle [matrix-basis-cos-angle
((Matrix Number) (Matrix Number) -> Number)] ((Matrix Number) (Matrix Number) -> Number)]
[matrix-basis-angle [matrix-basis-angle
@ -167,6 +182,10 @@
;; matrix-constructors.rkt ;; matrix-constructors.rkt
vandermonde-matrix vandermonde-matrix
;; matrix-basic.rkt ;; matrix-basic.rkt
matrix-1norm
matrix-2norm
matrix-inf-norm
matrix-norm
matrix-dot matrix-dot
matrix-cos-angle matrix-cos-angle
matrix-angle matrix-angle
@ -179,6 +198,9 @@
;; matrix-subspace.rkt ;; matrix-subspace.rkt
matrix-col-space matrix-col-space
;; matrix-operator-norm.rkt ;; matrix-operator-norm.rkt
matrix-op-1norm
matrix-op-2norm
matrix-op-inf-norm
matrix-basis-cos-angle matrix-basis-cos-angle
matrix-basis-angle matrix-basis-angle
;; matrix-qr.rkt ;; matrix-qr.rkt

173
collects/math/private/matrix/matrix-basic.rkt
View File

@ -130,29 +130,37 @@
(define i (unsafe-vector-ref js 0)) (define i (unsafe-vector-ref js 0))
(proc ((inst vector Index) i i)))))) (proc ((inst vector Index) i i))))))
(: matrix-upper-triangle (All (A) ((Matrix A) -> (Matrix (U A 0))))) (: matrix-upper-triangle (All (A) (case-> ((Matrix A) -> (Matrix (U A 0)))
(define (matrix-upper-triangle M) ((Matrix A) A -> (Matrix A)))))
(define-values (m n) (matrix-shape M)) (define matrix-upper-triangle
(define proc (unsafe-array-proc M)) (case-lambda
(array-default-strict [(M) (matrix-upper-triangle M 0)]
(unsafe-build-array [(M zero)
((inst vector Index) m n) (define-values (m n) (matrix-shape M))
(λ: ([ij : Indexes]) (define proc (unsafe-array-proc M))
(define i (unsafe-vector-ref ij 0)) (array-default-strict
(define j (unsafe-vector-ref ij 1)) (unsafe-build-array
(if (i . fx<= . j) (proc ij) 0))))) ((inst vector Index) m n)
(λ: ([ij : Indexes])
(define i (unsafe-vector-ref ij 0))
(define j (unsafe-vector-ref ij 1))
(if (i . fx<= . j) (proc ij) zero))))]))
(: matrix-lower-triangle (All (A) ((Matrix A) -> (Matrix (U A 0))))) (: matrix-lower-triangle (All (A) (case-> ((Matrix A) -> (Matrix (U A 0)))
(define (matrix-lower-triangle M) ((Matrix A) A -> (Matrix A)))))
(define-values (m n) (matrix-shape M)) (define matrix-lower-triangle
(define proc (unsafe-array-proc M)) (case-lambda
(array-default-strict [(M) (matrix-lower-triangle M 0)]
(unsafe-build-array [(M zero)
((inst vector Index) m n) (define-values (m n) (matrix-shape M))
(λ: ([ij : Indexes]) (define proc (unsafe-array-proc M))
(define i (unsafe-vector-ref ij 0)) (array-default-strict
(define j (unsafe-vector-ref ij 1)) (unsafe-build-array
(if (i . fx>= . j) (proc ij) 0))))) ((inst vector Index) m n)
(λ: ([ij : Indexes])
(define i (unsafe-vector-ref ij 0))
(define j (unsafe-vector-ref ij 1))
(if (i . fx>= . j) (proc ij) zero))))]))
;; =================================================================================================== ;; ===================================================================================================
;; Embiggenment (this is a perfectly cromulent word) ;; Embiggenment (this is a perfectly cromulent word)
@ -184,42 +192,63 @@
;; =================================================================================================== ;; ===================================================================================================
;; Inner product space (entrywise norm) ;; Inner product space (entrywise norm)
(: matrix-1norm ((Matrix Number) -> Nonnegative-Real)) (: nonstupid-magnitude (case-> (Flonum -> Nonnegative-Flonum)
(define (matrix-1norm a) (Real -> Nonnegative-Real)
(parameterize ([array-strictness #f]) (Float-Complex -> Nonnegative-Flonum)
(array-all-sum (array-magnitude a)))) (Number -> Nonnegative-Real)))
(define (nonstupid-magnitude x)
(if (real? x) (abs x) (magnitude x)))
(: matrix-2norm ((Matrix Number) -> Nonnegative-Real)) (: matrix-1norm (case-> ((Matrix Flonum) -> Nonnegative-Flonum)
(define (matrix-2norm a) ((Matrix Real) -> Nonnegative-Real)
((Matrix Float-Complex) -> Nonnegative-Flonum)
((Matrix Number) -> Nonnegative-Real)))
(define (matrix-1norm M)
(parameterize ([array-strictness #f]) (parameterize ([array-strictness #f])
(let ([a (array-strict (array-magnitude a))]) (array-all-sum (inline-array-map nonstupid-magnitude M))))
(: matrix-2norm (case-> ((Matrix Flonum) -> Nonnegative-Flonum)
((Matrix Real) -> Nonnegative-Real)
((Matrix Float-Complex) -> Nonnegative-Flonum)
((Matrix Number) -> Nonnegative-Real)))
(define (matrix-2norm M)
(parameterize ([array-strictness #f])
(let ([M (array-strict (inline-array-map nonstupid-magnitude M))])
;; Compute this divided by the maximum to avoid underflow and overflow ;; Compute this divided by the maximum to avoid underflow and overflow
(define mx (array-all-max a)) (define mx (array-all-max M))
(cond [(and (rational? mx) (positive? mx)) (cond [(and (rational? mx) (positive? mx))
(* mx (sqrt (array-all-sum (* mx (sqrt (array-all-sum (inline-array-map (λ (x) (sqr (/ x mx))) M))))]
(inline-array-map (λ: ([x : Nonnegative-Real]) (sqr (/ x mx))) a))))]
[else mx])))) [else mx]))))
(: matrix-inf-norm ((Matrix Number) -> Nonnegative-Real)) (: matrix-inf-norm (case-> ((Matrix Flonum) -> Nonnegative-Flonum)
(define (matrix-inf-norm a) ((Matrix Real) -> Nonnegative-Real)
((Matrix Float-Complex) -> Nonnegative-Flonum)
((Matrix Number) -> Nonnegative-Real)))
(define (matrix-inf-norm M)
(parameterize ([array-strictness #f]) (parameterize ([array-strictness #f])
(array-all-max (array-magnitude a)))) (array-all-max (inline-array-map nonstupid-magnitude M))))
(: matrix-p-norm ((Matrix Number) Positive-Real -> Nonnegative-Real)) (: matrix-p-norm (case-> ((Matrix Flonum) Positive-Real -> Nonnegative-Flonum)
(define (matrix-p-norm a p) ((Matrix Real) Positive-Real -> Nonnegative-Real)
((Matrix Float-Complex) Positive-Real -> Nonnegative-Flonum)
((Matrix Number) Positive-Real -> Nonnegative-Real)))
(define (matrix-p-norm M p)
(parameterize ([array-strictness #f]) (parameterize ([array-strictness #f])
(let ([a (array-strict (array-magnitude a))]) (let ([M (array-strict (inline-array-map nonstupid-magnitude M))])
;; Compute this divided by the maximum to avoid underflow and overflow ;; Compute this divided by the maximum to avoid underflow and overflow
(define mx (array-all-max a)) (define mx (array-all-max M))
(cond [(and (rational? mx) (positive? mx)) (cond [(and (rational? mx) (positive? mx))
(assert (* mx (expt (array-all-sum (inline-array-map (λ (x) (expt (abs (/ x mx)) p)) M))
(* mx (expt (array-all-sum (/ p)))]
(inline-array-map (λ: ([x : Nonnegative-Real]) (expt (/ x mx) p)) a))
(/ p)))
(make-predicate Nonnegative-Real))]
[else mx])))) [else mx]))))
(: matrix-norm (case-> ((Matrix Number) -> Nonnegative-Real) (: matrix-norm (case-> ((Matrix Flonum) -> Nonnegative-Flonum)
((Matrix Flonum) Real -> Nonnegative-Flonum)
((Matrix Real) -> Nonnegative-Real)
((Matrix Real) Real -> Nonnegative-Real)
((Matrix Float-Complex) -> Nonnegative-Flonum)
((Matrix Float-Complex) Real -> Nonnegative-Flonum)
((Matrix Number) -> Nonnegative-Real)
((Matrix Number) Real -> Nonnegative-Real))) ((Matrix Number) Real -> Nonnegative-Real)))
;; Computes the p norm of a matrix ;; Computes the p norm of a matrix
(define (matrix-norm a [p 2]) (define (matrix-norm a [p 2])
@ -230,8 +259,12 @@
[(p . > . 1) (matrix-p-norm a p)] [(p . > . 1) (matrix-p-norm a p)]
[else (raise-argument-error 'matrix-norm "Real >= 1" 1 a p)])) [else (raise-argument-error 'matrix-norm "Real >= 1" 1 a p)]))
(: matrix-dot (case-> ((Matrix Real) -> Nonnegative-Real) (: matrix-dot (case-> ((Matrix Flonum) -> Nonnegative-Flonum)
((Matrix Flonum) (Matrix Flonum) -> Flonum)
((Matrix Real) -> Nonnegative-Real)
((Matrix Real) (Matrix Real) -> Real) ((Matrix Real) (Matrix Real) -> Real)
((Matrix Float-Complex) -> Nonnegative-Flonum)
((Matrix Float-Complex) (Matrix Float-Complex) -> Float-Complex)
((Matrix Number) -> Nonnegative-Real) ((Matrix Number) -> Nonnegative-Real)
((Matrix Number) (Matrix Number) -> Number))) ((Matrix Number) (Matrix Number) -> Number)))
;; Computes the Frobenius inner product of a matrix with itself or of two matrices ;; Computes the Frobenius inner product of a matrix with itself or of two matrices
@ -256,20 +289,30 @@
(λ: ([js : Indexes]) (λ: ([js : Indexes])
(* (aproc js) (conjugate (bproc js)))))))])) (* (aproc js) (conjugate (bproc js)))))))]))
(: matrix-cos-angle (case-> ((Matrix Real) (Matrix Real) -> Real) (: matrix-cos-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))) ((Matrix Number) (Matrix Number) -> Number)))
(define (matrix-cos-angle M N) (define (matrix-cos-angle M N)
(/ (matrix-dot M N) (* (matrix-2norm M) (matrix-2norm N)))) (/ (matrix-dot M N) (* (matrix-2norm M) (matrix-2norm N))))
(: matrix-angle (case-> ((Matrix Real) (Matrix Real) -> Real) (: 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))) ((Matrix Number) (Matrix Number) -> Number)))
(define (matrix-angle M N) (define (matrix-angle M N)
(acos (matrix-cos-angle M N))) (acos (matrix-cos-angle M N)))
(: matrix-normalize (: matrix-normalize
(All (A) (case-> ((Matrix Real) -> (Matrix Real)) (All (A) (case-> ((Matrix Flonum) -> (Matrix Flonum))
((Matrix Flonum) Real -> (Matrix Flonum))
((Matrix Flonum) Real (-> A) -> (U A (Matrix Flonum)))
((Matrix Real) -> (Matrix Real))
((Matrix Real) Real -> (Matrix Real)) ((Matrix Real) Real -> (Matrix Real))
((Matrix Real) Real (-> A) -> (U A (Matrix Real))) ((Matrix Real) Real (-> A) -> (U A (Matrix Real)))
((Matrix Float-Complex) -> (Matrix Float-Complex))
((Matrix Float-Complex) Real -> (Matrix Float-Complex))
((Matrix Float-Complex) Real (-> A) -> (U A (Matrix Float-Complex)))
((Matrix Number) -> (Matrix Number)) ((Matrix Number) -> (Matrix Number))
((Matrix Number) Real -> (Matrix Number)) ((Matrix Number) Real -> (Matrix Number))
((Matrix Number) Real (-> A) -> (U A (Matrix Number)))))) ((Matrix Number) Real (-> A) -> (U A (Matrix Number))))))
@ -291,19 +334,25 @@
(define (matrix-transpose a) (define (matrix-transpose a)
(array-axis-swap (ensure-matrix 'matrix-transpose a) 0 1)) (array-axis-swap (ensure-matrix 'matrix-transpose a) 0 1))
(: matrix-conjugate (case-> ((Matrix Real) -> (Matrix Real)) (: matrix-conjugate (case-> ((Matrix Flonum) -> (Matrix Flonum))
((Matrix Real) -> (Matrix Real))
((Matrix Float-Complex) -> (Matrix Float-Complex))
((Matrix Number) -> (Matrix Number)))) ((Matrix Number) -> (Matrix Number))))
(define (matrix-conjugate a) (define (matrix-conjugate a)
(array-conjugate (ensure-matrix 'matrix-conjugate a))) (array-conjugate (ensure-matrix 'matrix-conjugate a)))
(: matrix-hermitian (case-> ((Matrix Real) -> (Matrix Real)) (: matrix-hermitian (case-> ((Matrix Flonum) -> (Matrix Flonum))
((Matrix Real) -> (Matrix Real))
((Matrix Float-Complex) -> (Matrix Float-Complex))
((Matrix Number) -> (Matrix Number)))) ((Matrix Number) -> (Matrix Number))))
(define (matrix-hermitian a) (define (matrix-hermitian a)
(array-default-strict (array-default-strict
(parameterize ([array-strictness #f]) (parameterize ([array-strictness #f])
(array-axis-swap (array-conjugate (ensure-matrix 'matrix-hermitian a)) 0 1)))) (array-axis-swap (array-conjugate (ensure-matrix 'matrix-hermitian a)) 0 1))))
(: matrix-trace (case-> ((Matrix Real) -> Real) (: matrix-trace (case-> ((Matrix Flonum) -> Flonum)
((Matrix Real) -> Real)
((Matrix Float-Complex) -> Float-Complex)
((Matrix Number) -> Number))) ((Matrix Number) -> Number)))
(define (matrix-trace a) (define (matrix-trace a)
(cond [(square-matrix? a) (cond [(square-matrix? a)
@ -349,15 +398,23 @@
[else (fail)])]))] [else (fail)])]))]
[else (fail)])])) [else (fail)])]))
(: make-matrix-normalize (Real -> (case-> ((Matrix Real) -> (U #f (Matrix Real))) (: make-matrix-normalize (Real -> (case-> ((Matrix Flonum) -> (U #f (Matrix Flonum)))
((Matrix Real) -> (U #f (Matrix Real)))
((Matrix Float-Complex) -> (U #f (Matrix Float-Complex)))
((Matrix Number) -> (U #f (Matrix Number)))))) ((Matrix Number) -> (U #f (Matrix Number))))))
(define ((make-matrix-normalize p) M) (define ((make-matrix-normalize p) M)
(matrix-normalize M p (λ () #f))) (matrix-normalize M p (λ () #f)))
(: matrix-normalize-rows (: matrix-normalize-rows
(All (A) (case-> ((Matrix Real) -> (Matrix Real)) (All (A) (case-> ((Matrix Flonum) -> (Matrix Flonum))
((Matrix Flonum) Real -> (Matrix Flonum))
((Matrix Flonum) Real (-> A) -> (U A (Matrix Flonum)))
((Matrix Real) -> (Matrix Real))
((Matrix Real) Real -> (Matrix Real)) ((Matrix Real) Real -> (Matrix Real))
((Matrix Real) Real (-> A) -> (U A (Matrix Real))) ((Matrix Real) Real (-> A) -> (U A (Matrix Real)))
((Matrix Float-Complex) -> (Matrix Float-Complex))
((Matrix Float-Complex) Real -> (Matrix Float-Complex))
((Matrix Float-Complex) Real (-> A) -> (U A (Matrix Float-Complex)))
((Matrix Number) -> (Matrix Number)) ((Matrix Number) -> (Matrix Number))
((Matrix Number) Real -> (Matrix Number)) ((Matrix Number) Real -> (Matrix Number))
((Matrix Number) Real (-> A) -> (U A (Matrix Number)))))) ((Matrix Number) Real (-> A) -> (U A (Matrix Number))))))
@ -371,9 +428,15 @@
(matrix-map-rows (make-matrix-normalize p) M fail)])) (matrix-map-rows (make-matrix-normalize p) M fail)]))
(: matrix-normalize-cols (: matrix-normalize-cols
(All (A) (case-> ((Matrix Real) -> (Matrix Real)) (All (A) (case-> ((Matrix Flonum) -> (Matrix Flonum))
((Matrix Flonum) Real -> (Matrix Flonum))
((Matrix Flonum) Real (-> A) -> (U A (Matrix Flonum)))
((Matrix Real) -> (Matrix Real))
((Matrix Real) Real -> (Matrix Real)) ((Matrix Real) Real -> (Matrix Real))
((Matrix Real) Real (-> A) -> (U A (Matrix Real))) ((Matrix Real) Real (-> A) -> (U A (Matrix Real)))
((Matrix Float-Complex) -> (Matrix Float-Complex))
((Matrix Float-Complex) Real -> (Matrix Float-Complex))
((Matrix Float-Complex) Real (-> A) -> (U A (Matrix Float-Complex)))
((Matrix Number) -> (Matrix Number)) ((Matrix Number) -> (Matrix Number))
((Matrix Number) Real -> (Matrix Number)) ((Matrix Number) Real -> (Matrix Number))
((Matrix Number) Real (-> A) -> (U A (Matrix Number)))))) ((Matrix Number) Real (-> A) -> (U A (Matrix Number))))))

49
collects/math/private/matrix/matrix-constructors.rkt
View File

@ -1,8 +1,10 @@
#lang typed/racket/base #lang typed/racket/base
(require racket/fixnum (require racket/fixnum
racket/flonum
racket/list racket/list
racket/vector racket/vector
math/base
"matrix-types.rkt" "matrix-types.rkt"
"../unsafe.rkt" "../unsafe.rkt"
"../array/array-struct.rkt" "../array/array-struct.rkt"
@ -22,8 +24,14 @@
;; =================================================================================================== ;; ===================================================================================================
;; Basic constructors ;; Basic constructors
(: identity-matrix (Integer -> (Matrix (U 0 1)))) (: identity-matrix (All (A) (case-> (Integer -> (Matrix (U 1 0)))
(define (identity-matrix m) (diagonal-array 2 m 1 0)) (Integer A -> (Matrix (U A 0)))
(Integer A A -> (Matrix A)))))
(define identity-matrix
(case-lambda
[(m) (diagonal-array 2 m 1 0)]
[(m one) (diagonal-array 2 m one 0)]
[(m one zero) (diagonal-array 2 m one zero)]))
(: make-matrix (All (A) (Integer Integer A -> (Matrix A)))) (: make-matrix (All (A) (Integer Integer A -> (Matrix A))))
(define (make-matrix m n x) (define (make-matrix m n x)
@ -60,9 +68,12 @@
(cond [(= i (unsafe-vector-ref js 1)) (unsafe-vector-ref vs i)] (cond [(= i (unsafe-vector-ref js 1)) (unsafe-vector-ref vs i)]
[else zero])))])) [else zero])))]))
(: diagonal-matrix (All (A) ((Listof A) -> (Matrix (U A 0))))) (: diagonal-matrix (All (A) (case-> ((Listof A) -> (Matrix (U A 0)))
(define (diagonal-matrix xs) ((Listof A) A -> (Matrix A)))))
(diagonal-matrix/zero xs 0)) (define diagonal-matrix
(case-lambda
[(xs) (diagonal-matrix/zero xs 0)]
[(xs zero) (diagonal-matrix/zero xs zero)]))
;; =================================================================================================== ;; ===================================================================================================
;; Block diagonal matrices ;; Block diagonal matrices
@ -129,21 +140,29 @@
[else [else
(block-diagonal-matrix/zero* as zero)]))) (block-diagonal-matrix/zero* as zero)])))
(: block-diagonal-matrix (All (A) ((Listof (Matrix A)) -> (Matrix (U A 0))))) (: block-diagonal-matrix (All (A) (case-> ((Listof (Matrix A)) -> (Matrix (U A 0)))
(define (block-diagonal-matrix as) ((Listof (Matrix A)) A -> (Matrix A)))))
(block-diagonal-matrix/zero as 0)) (define block-diagonal-matrix
(case-lambda
[(as) (block-diagonal-matrix/zero as 0)]
[(as zero) (block-diagonal-matrix/zero as zero)]))
;; =================================================================================================== ;; ===================================================================================================
;; Special matrices ;; Special matrices
(: expt-hack (case-> (Real Integer -> Real) (: sane-expt (case-> (Flonum Index -> Flonum)
(Number Integer -> Number))) (Real Index -> Real)
;; Stop using this when TR correctly derives expt : Real Integer -> Real (Float-Complex Index -> Float-Complex)
(define (expt-hack x n) (Number Index -> Number)))
(cond [(real? x) (assert (expt x n) real?)] (define (sane-expt x n)
(cond [(flonum? x) (flexpt x (real->double-flonum n))]
[(real? x) (real-part (expt x n))] ; remove `real-part' when expt : Real Index -> Real
[(float-complex? x) (number->float-complex (expt x n))]
[else (expt x n)])) [else (expt x n)]))
(: vandermonde-matrix (case-> ((Listof Real) Integer -> (Matrix Real)) (: vandermonde-matrix (case-> ((Listof Flonum) Integer -> (Matrix Flonum))
((Listof Real) Integer -> (Matrix Real))
((Listof Float-Complex) Integer -> (Matrix Float-Complex))
((Listof Number) Integer -> (Matrix Number)))) ((Listof Number) Integer -> (Matrix Number))))
(define (vandermonde-matrix xs n) (define (vandermonde-matrix xs n)
(cond [(empty? xs) (cond [(empty? xs)
@ -151,4 +170,4 @@
[(or (not (index? n)) (zero? n)) [(or (not (index? n)) (zero? n))
(raise-argument-error 'vandermonde-matrix "Positive-Index" 1 xs n)] (raise-argument-error 'vandermonde-matrix "Positive-Index" 1 xs n)]
[else [else
(array-axis-expand (list->array xs) 1 n expt-hack)])) (array-axis-expand (list->array xs) 1 n sane-expt)]))

22
collects/math/private/matrix/matrix-expt.rkt
View File

@ -3,11 +3,16 @@
(require "matrix-types.rkt" (require "matrix-types.rkt"
"matrix-constructors.rkt" "matrix-constructors.rkt"
"matrix-arithmetic.rkt" "matrix-arithmetic.rkt"
"utils.rkt") "utils.rkt"
"../array/array-struct.rkt"
"../array/utils.rkt"
"../unsafe.rkt")
(provide matrix-expt) (provide matrix-expt)
(: matrix-expt/ns (case-> ((Matrix Real) Positive-Integer -> (Matrix Real)) (: matrix-expt/ns (case-> ((Matrix Flonum) Positive-Integer -> (Matrix Flonum))
((Matrix Real) Positive-Integer -> (Matrix Real))
((Matrix Float-Complex) Positive-Integer -> (Matrix Float-Complex))
((Matrix Number) Positive-Integer -> (Matrix Number)))) ((Matrix Number) Positive-Integer -> (Matrix Number))))
(define (matrix-expt/ns a n) (define (matrix-expt/ns a n)
(define n/2 (quotient n 2)) (define n/2 (quotient n 2))
@ -22,10 +27,19 @@
;; m = n - 1 ;; m = n - 1
(matrix* a (matrix-expt/ns a m)))))) (matrix* a (matrix-expt/ns a m))))))
(: matrix-expt (case-> ((Matrix Real) Integer -> (Matrix Real)) (: matrix-expt (case-> ((Matrix Flonum) Integer -> (Matrix Flonum))
((Matrix Real) Integer -> (Matrix Real))
((Matrix Float-Complex) Integer -> (Matrix Float-Complex))
((Matrix Number) Integer -> (Matrix Number)))) ((Matrix Number) Integer -> (Matrix Number))))
(define (matrix-expt a n) (define (matrix-expt a n)
(cond [(not (square-matrix? a)) (raise-argument-error 'matrix-expt "square-matrix?" 0 a n)] (cond [(not (square-matrix? a)) (raise-argument-error 'matrix-expt "square-matrix?" 0 a n)]
[(negative? n) (raise-argument-error 'matrix-expt "Natural" 1 a n)] [(negative? n) (raise-argument-error 'matrix-expt "Natural" 1 a n)]
[(zero? n) (identity-matrix (square-matrix-size a))] [(zero? n)
(define proc (unsafe-array-proc a))
(array-default-strict
(unsafe-build-array
(array-shape a)
(λ: ([ij : Indexes])
(define x (proc ij))
(if (= (unsafe-vector-ref ij 0) (unsafe-vector-ref ij 1)) (one* x) (zero* x)))))]
[else (call/ns (λ () (matrix-expt/ns a n)))])) [else (call/ns (λ () (matrix-expt/ns a n)))]))

24
collects/math/private/matrix/matrix-gauss-elim.rkt
View File

@ -16,10 +16,18 @@
(define-type Pivoting (U 'first 'partial)) (define-type Pivoting (U 'first 'partial))
(: matrix-gauss-elim (: matrix-gauss-elim
(case-> ((Matrix Real) -> (Values (Matrix Real) (Listof Index))) (case-> ((Matrix Flonum) -> (Values (Matrix Flonum) (Listof Index)))
((Matrix Flonum) Any -> (Values (Matrix Flonum) (Listof Index)))
((Matrix Flonum) Any Any -> (Values (Matrix Flonum) (Listof Index)))
((Matrix Flonum) Any Any Pivoting -> (Values (Matrix Flonum) (Listof Index)))
((Matrix Real) -> (Values (Matrix Real) (Listof Index)))
((Matrix Real) Any -> (Values (Matrix Real) (Listof Index))) ((Matrix Real) Any -> (Values (Matrix Real) (Listof Index)))
((Matrix Real) Any Any -> (Values (Matrix Real) (Listof Index))) ((Matrix Real) Any Any -> (Values (Matrix Real) (Listof Index)))
((Matrix Real) Any Any Pivoting -> (Values (Matrix Real) (Listof Index))) ((Matrix Real) Any Any Pivoting -> (Values (Matrix Real) (Listof Index)))
((Matrix Float-Complex) -> (Values (Matrix Float-Complex) (Listof Index)))
((Matrix Float-Complex) Any -> (Values (Matrix Float-Complex) (Listof Index)))
((Matrix Float-Complex) Any Any -> (Values (Matrix Float-Complex) (Listof Index)))
((Matrix Float-Complex) Any Any Pivoting -> (Values (Matrix Float-Complex) (Listof Index)))
((Matrix Number) -> (Values (Matrix Number) (Listof Index))) ((Matrix Number) -> (Values (Matrix Number) (Listof Index)))
((Matrix Number) Any -> (Values (Matrix Number) (Listof Index))) ((Matrix Number) Any -> (Values (Matrix Number) (Listof Index)))
((Matrix Number) Any Any -> (Values (Matrix Number) (Listof Index))) ((Matrix Number) Any Any -> (Values (Matrix Number) (Listof Index)))
@ -52,17 +60,25 @@
(vector-swap! rows i p) (vector-swap! rows i p)
;; Possibly unitize the new current row ;; Possibly unitize the new current row
(let ([pivot (if unitize-pivot? (let ([pivot (if unitize-pivot?
(begin (vector-scale! (unsafe-vector-ref rows i) (/ pivot)) (begin (vector-scale! (unsafe-vector-ref rows i) (/ 1 pivot))
1) (/ pivot pivot))
pivot)]) pivot)])
(elim-rows! rows m i j pivot (if jordan? 0 (fx+ i 1))) (elim-rows! rows m i j pivot (if jordan? 0 (fx+ i 1)))
(loop (fx+ i 1) (fx+ j 1) without-pivot))])]))) (loop (fx+ i 1) (fx+ j 1) without-pivot))])])))
(: matrix-row-echelon (: matrix-row-echelon
(case-> ((Matrix Real) -> (Matrix Real)) (case-> ((Matrix Flonum) -> (Matrix Flonum))
((Matrix Flonum) Any -> (Matrix Flonum))
((Matrix Flonum) Any Any -> (Matrix Flonum))
((Matrix Flonum) Any Any Pivoting -> (Matrix Flonum))
((Matrix Real) -> (Matrix Real))
((Matrix Real) Any -> (Matrix Real)) ((Matrix Real) Any -> (Matrix Real))
((Matrix Real) Any Any -> (Matrix Real)) ((Matrix Real) Any Any -> (Matrix Real))
((Matrix Real) Any Any Pivoting -> (Matrix Real)) ((Matrix Real) Any Any Pivoting -> (Matrix Real))
((Matrix Float-Complex) -> (Matrix Float-Complex))
((Matrix Float-Complex) Any -> (Matrix Float-Complex))
((Matrix Float-Complex) Any Any -> (Matrix Float-Complex))
((Matrix Float-Complex) Any Any Pivoting -> (Matrix Float-Complex))
((Matrix Number) -> (Matrix Number)) ((Matrix Number) -> (Matrix Number))
((Matrix Number) Any -> (Matrix Number)) ((Matrix Number) Any -> (Matrix Number))
((Matrix Number) Any Any -> (Matrix Number)) ((Matrix Number) Any Any -> (Matrix Number))

40
collects/math/private/matrix/matrix-gram-schmidt.rkt
View File

@ -16,7 +16,9 @@
(provide matrix-gram-schmidt (provide matrix-gram-schmidt
matrix-basis-extension) matrix-basis-extension)
(: find-nonzero-vector (case-> ((Vectorof (Vectorof Real)) -> (U #f Index)) (: find-nonzero-vector (case-> ((Vectorof (Vectorof Flonum)) -> (U #f Index))
((Vectorof (Vectorof Real)) -> (U #f Index))
((Vectorof (Vectorof Float-Complex)) -> (U #f Index))
((Vectorof (Vectorof Number)) -> (U #f Index)))) ((Vectorof (Vectorof Number)) -> (U #f Index))))
(define (find-nonzero-vector vss) (define (find-nonzero-vector vss)
(define n (vector-length vss)) (define n (vector-length vss))
@ -28,7 +30,10 @@
[else #f]))])) [else #f]))]))
(: subtract-projections! (: subtract-projections!
(case-> ((Vectorof (Vectorof Real)) Nonnegative-Fixnum Index (Vectorof Real) -> Void) (case-> ((Vectorof (Vectorof Flonum)) Nonnegative-Fixnum Index (Vectorof Flonum) -> Void)
((Vectorof (Vectorof Real)) Nonnegative-Fixnum Index (Vectorof Real) -> Void)
((Vectorof (Vectorof Float-Complex)) Nonnegative-Fixnum Index (Vectorof Float-Complex)
-> Void)
((Vectorof (Vectorof Number)) Nonnegative-Fixnum Index (Vectorof Number) -> Void))) ((Vectorof (Vectorof Number)) Nonnegative-Fixnum Index (Vectorof Number) -> Void)))
(define (subtract-projections! rows i m row) (define (subtract-projections! rows i m row)
(let loop ([#{i : Nonnegative-Fixnum} i]) (let loop ([#{i : Nonnegative-Fixnum} i])
@ -36,7 +41,9 @@
(vector-sub-proj! (unsafe-vector-ref rows i) row #f) (vector-sub-proj! (unsafe-vector-ref rows i) row #f)
(loop (fx+ i 1))))) (loop (fx+ i 1)))))
(: matrix-gram-schmidt/ns (case-> ((Matrix Real) Any Integer -> (Array Real)) (: matrix-gram-schmidt/ns (case-> ((Matrix Flonum) Any Integer -> (Array Flonum))
((Matrix Real) Any Integer -> (Array Real))
((Matrix Float-Complex) Any Integer -> (Array Float-Complex))
((Matrix Number) Any Integer -> (Array Number)))) ((Matrix Number) Any Integer -> (Array Number))))
;; Performs Gram-Schmidt orthogonalization on M, assuming the rows before `start' are already ;; Performs Gram-Schmidt orthogonalization on M, assuming the rows before `start' are already
;; orthogonal ;; orthogonal
@ -60,31 +67,46 @@
[else [else
(matrix-transpose (vector*->matrix (list->vector (reverse bs))))]))] (matrix-transpose (vector*->matrix (list->vector (reverse bs))))]))]
[else [else
(make-array (vector (matrix-num-rows M) 0) 0)])) (make-array (vector (matrix-num-rows M) 0)
;; Value won't be in the matrix, but this satisfies TR:
(zero* (unsafe-vector2d-ref rows 0 0)))]))
(: matrix-gram-schmidt (case-> ((Matrix Real) -> (Array Real)) (: matrix-gram-schmidt (case-> ((Matrix Flonum) -> (Array Flonum))
((Matrix Flonum) Any -> (Array Flonum))
((Matrix Flonum) Any Integer -> (Array Flonum))
((Matrix Real) -> (Array Real))
((Matrix Real) Any -> (Array Real)) ((Matrix Real) Any -> (Array Real))
((Matrix Real) Any Integer -> (Array Real)) ((Matrix Real) Any Integer -> (Array Real))
((Matrix Float-Complex) -> (Array Float-Complex))
((Matrix Float-Complex) Any -> (Array Float-Complex))
((Matrix Float-Complex) Any Integer -> (Array Float-Complex))
((Matrix Number) -> (Array Number)) ((Matrix Number) -> (Array Number))
((Matrix Number) Any -> (Array Number)) ((Matrix Number) Any -> (Array Number))
((Matrix Number) Any Integer -> (Array Number)))) ((Matrix Number) Any Integer -> (Array Number))))
(define (matrix-gram-schmidt M [normalize? #f] [start 0]) (define (matrix-gram-schmidt M [normalize? #f] [start 0])
(call/ns (λ () (matrix-gram-schmidt/ns M normalize? start)))) (call/ns (λ () (matrix-gram-schmidt/ns M normalize? start))))
(: matrix-basis-extension/ns (case-> ((Matrix Real) -> (Array Real)) (: matrix-basis-extension/ns (case-> ((Matrix Flonum) -> (Array Flonum))
((Matrix Real) -> (Array Real))
((Matrix Float-Complex) -> (Array Float-Complex))
((Matrix Number) -> (Array Number)))) ((Matrix Number) -> (Array Number))))
(define (matrix-basis-extension/ns B) (define (matrix-basis-extension/ns B)
(define-values (m n) (matrix-shape B)) (define-values (m n) (matrix-shape B))
(define x00 (matrix-ref B 0 0))
(define zero (zero* x00))
(define one (one* x00))
(cond [(n . < . m) (cond [(n . < . m)
(define S (matrix-gram-schmidt (matrix-augment (list B (identity-matrix m))) #f n)) (define S (matrix-gram-schmidt (matrix-augment (list B (identity-matrix m one zero))) #f n))
(define R (submatrix S (::) (:: n #f))) (define R (submatrix S (::) (:: n #f)))
(matrix-augment (take (sort/key (matrix-cols R) > matrix-norm) (- m n)))] (matrix-augment (take (sort/key (matrix-cols R) > matrix-norm) (- m n)))]
[(n . = . m) [(n . = . m)
(make-array (vector m 0) 0)] (make-array (vector m 0) zero)]
[else [else
(raise-argument-error 'matrix-extend-row-basis "matrix? with width < height" B)])) (raise-argument-error 'matrix-extend-row-basis "matrix? with width < height" B)]))
(: matrix-basis-extension (case-> ((Matrix Real) -> (Array Real)) (: matrix-basis-extension (case-> ((Matrix Flonum) -> (Array Flonum))
((Matrix Real) -> (Array Real))
((Matrix Float-Complex) -> (Array Float-Complex))
((Matrix Number) -> (Array Number)))) ((Matrix Number) -> (Array Number))))
(define (matrix-basis-extension B) (define (matrix-basis-extension B)
(call/ns (λ () (matrix-basis-extension/ns B)))) (call/ns (λ () (matrix-basis-extension/ns B))))

18
collects/math/private/matrix/matrix-lu.rkt
View File

@ -8,15 +8,22 @@
"../unsafe.rkt" "../unsafe.rkt"
"../vector/vector-mutate.rkt" "../vector/vector-mutate.rkt"
"../array/mutable-array.rkt" "../array/mutable-array.rkt"
"../array/array-struct.rkt") "../array/array-struct.rkt"
"../array/array-pointwise.rkt")
(provide matrix-lu) (provide matrix-lu)
;; An LU factorization exists iff Gaussian elimination can be done without row swaps. ;; An LU factorization exists iff Gaussian elimination can be done without row swaps.
(: matrix-lu (: matrix-lu
(All (A) (case-> ((Matrix Real) -> (Values (Matrix Real) (Matrix Real))) (All (A) (case-> ((Matrix Flonum) -> (Values (Matrix Flonum) (Matrix Flonum)))
((Matrix Flonum) (-> A) -> (Values (U A (Matrix Flonum)) (Matrix Flonum)))
((Matrix Real) -> (Values (Matrix Real) (Matrix Real)))
((Matrix Real) (-> A) -> (Values (U A (Matrix Real)) (Matrix Real))) ((Matrix Real) (-> A) -> (Values (U A (Matrix Real)) (Matrix Real)))
((Matrix Float-Complex) -> (Values (Matrix Float-Complex)
(Matrix Float-Complex)))
((Matrix Float-Complex) (-> A) -> (Values (U A (Matrix Float-Complex))
(Matrix Float-Complex)))
((Matrix Number) -> (Values (Matrix Number) (Matrix Number))) ((Matrix Number) -> (Values (Matrix Number) (Matrix Number)))
((Matrix Number) (-> A) -> (Values (U A (Matrix Number)) (Matrix Number)))))) ((Matrix Number) (-> A) -> (Values (U A (Matrix Number)) (Matrix Number))))))
(define matrix-lu (define matrix-lu
@ -28,7 +35,7 @@
(define L (define L
(parameterize ([array-strictness #f]) (parameterize ([array-strictness #f])
;; Construct L in a weird way to prove to TR that it has the right type ;; Construct L in a weird way to prove to TR that it has the right type
(array->mutable-array (matrix-scale M (ann 0 Real))))) (array->mutable-array (inline-array-map zero* M))))
;; Going to fill in the lower triangle by banging values into `ys' ;; Going to fill in the lower triangle by banging values into `ys'
(define ys (mutable-array-data L)) (define ys (mutable-array-data L))
(let loop ([#{i : Nonnegative-Fixnum} 0]) (let loop ([#{i : Nonnegative-Fixnum} 0])
@ -51,12 +58,13 @@
;; Add row i, scaled ;; Add row i, scaled
(vector-scaled-add! (unsafe-vector-ref rows l) (vector-scaled-add! (unsafe-vector-ref rows l)
(unsafe-vector-ref rows i) (unsafe-vector-ref rows i)
(- y_li))) (* -1 y_li)))
(l-loop (fx+ l 1))] (l-loop (fx+ l 1))]
[else [else
(loop (fx+ i 1))]))])] (loop (fx+ i 1))]))])]
[else [else
;; L's lower triangle has been filled; now fill the diagonal with 1s ;; L's lower triangle has been filled; now fill the diagonal with 1s
(for: ([i : Integer (in-range 0 m)]) (for: ([i : Integer (in-range 0 m)])
(vector-set! ys (+ (* i m) i) 1)) (define j (+ (* i m) i))
(vector-set! ys j (one* (vector-ref ys j))))
(values L (vector*->matrix rows))]))])) (values L (vector*->matrix rows))]))]))

23
collects/math/private/matrix/matrix-operator-norm.rkt
View File

@ -42,33 +42,46 @@ See "How to Measure Errors" in the LAPACK manual for more details:
matrix-orthonormal? matrix-orthonormal?
) )
(: matrix-op-1norm ((Matrix Number) -> Nonnegative-Real)) (: matrix-op-1norm (case-> ((Matrix Flonum) -> Nonnegative-Flonum)
((Matrix Real) -> Nonnegative-Real)
((Matrix Float-Complex) -> Nonnegative-Flonum)
((Matrix Number) -> Nonnegative-Real)))
;; When M is a column matrix, this is equivalent to matrix-1norm ;; When M is a column matrix, this is equivalent to matrix-1norm
(define (matrix-op-1norm M) (define (matrix-op-1norm M)
(parameterize ([array-strictness #f]) (parameterize ([array-strictness #f])
(assert (apply max (map matrix-1norm (matrix-cols M))) nonnegative?))) (assert (apply max (map matrix-1norm (matrix-cols M))) nonnegative?)))
(: matrix-op-2norm ((Matrix Number) -> Nonnegative-Real)) (: matrix-op-2norm (case-> ((Matrix Flonum) -> Nonnegative-Flonum)
((Matrix Real) -> Nonnegative-Real)
((Matrix Float-Complex) -> Nonnegative-Flonum)
((Matrix Number) -> Nonnegative-Real)))
;; When M is a column matrix, this is equivalent to matrix-2norm ;; When M is a column matrix, this is equivalent to matrix-2norm
(define (matrix-op-2norm M) (define (matrix-op-2norm M)
;(matrix-max-singular-value M) ;(matrix-max-singular-value M)
;(sqrt (matrix-max-eigenvalue M)) ;(sqrt (matrix-max-eigenvalue M))
(error 'unimplemented)) (error 'unimplemented))
(: matrix-op-inf-norm ((Matrix Number) -> Nonnegative-Real)) (: matrix-op-inf-norm (case-> ((Matrix Flonum) -> Nonnegative-Flonum)
((Matrix Real) -> Nonnegative-Real)
((Matrix Float-Complex) -> Nonnegative-Flonum)
((Matrix Number) -> Nonnegative-Real)))
;; When M is a column matrix, this is equivalent to matrix-inf-norm ;; When M is a column matrix, this is equivalent to matrix-inf-norm
(define (matrix-op-inf-norm M) (define (matrix-op-inf-norm M)
(parameterize ([array-strictness #f]) (parameterize ([array-strictness #f])
(assert (apply max (map matrix-1norm (matrix-rows M))) nonnegative?))) (assert (apply max (map matrix-1norm (matrix-rows M))) nonnegative?)))
(: matrix-basis-cos-angle (case-> ((Matrix Real) (Matrix Real) -> Real) (: matrix-basis-cos-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))) ((Matrix Number) (Matrix Number) -> Number)))
;; Returns the angle between the two subspaces spanned by the two given sets of column vectors ;; Returns the angle between the two subspaces spanned by the two given sets of column vectors
(define (matrix-basis-cos-angle M R) (define (matrix-basis-cos-angle M R)
;(matrix-min-singular-value (matrix* (matrix-hermitian M) R)) ;(matrix-min-singular-value (matrix* (matrix-hermitian M) R))
(error 'unimplemented)) (error 'unimplemented))
(: matrix-basis-angle (case-> ((Matrix Real) (Matrix Real) -> Real) (: matrix-basis-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))) ((Matrix Number) (Matrix Number) -> Number)))
;; Returns the angle between the two subspaces spanned by the two given sets of column vectors ;; Returns the angle between the two subspaces spanned by the two given sets of column vectors
(define (matrix-basis-angle M R) (define (matrix-basis-angle M R)

23
collects/math/private/matrix/matrix-qr.rkt
View File

@ -5,6 +5,7 @@
"matrix-arithmetic.rkt" "matrix-arithmetic.rkt"
"matrix-constructors.rkt" "matrix-constructors.rkt"
"matrix-gram-schmidt.rkt" "matrix-gram-schmidt.rkt"
"utils.rkt"
"../array/array-transform.rkt" "../array/array-transform.rkt"
"../array/array-struct.rkt") "../array/array-struct.rkt")
@ -24,21 +25,33 @@ produces matrices for which `matrix-orthogonal?' returns #t with eps <= 10*epsil
independently of the matrix size. independently of the matrix size.
|# |#
(: matrix-qr/ns (case-> ((Matrix Real) Any -> (Values (Matrix Real) (Matrix Real))) (: matrix-qr/ns (case-> ((Matrix Flonum) Any -> (Values (Matrix Flonum) (Matrix Flonum)))
((Matrix Real) Any -> (Values (Matrix Real) (Matrix Real)))
((Matrix Float-Complex) Any -> (Values (Matrix Float-Complex)
(Matrix Float-Complex)))
((Matrix Number) Any -> (Values (Matrix Number) (Matrix Number))))) ((Matrix Number) Any -> (Values (Matrix Number) (Matrix Number)))))
(define (matrix-qr/ns M full?) (define (matrix-qr/ns M full?)
(define x00 (matrix-ref M 0 0))
(define zero (zero* x00))
(define one (one* x00))
(define B (matrix-gram-schmidt M #f)) (define B (matrix-gram-schmidt M #f))
(define Q (define Q
(matrix-gram-schmidt (matrix-gram-schmidt
(cond [(or (square-matrix? B) (and (matrix? B) (not full?))) B] (cond [(or (square-matrix? B) (and (matrix? B) (not full?))) B]
[(matrix? B) (array-append* (list B (matrix-basis-extension B)) 1)] [(matrix? B) (array-append* (list B (matrix-basis-extension B)) 1)]
[full? (identity-matrix (matrix-num-rows M))] [full? (identity-matrix (matrix-num-rows M) one zero)]
[else (matrix-col (identity-matrix (matrix-num-rows M)) 0)]) [else (matrix-col (identity-matrix (matrix-num-rows M) one zero) 0)])
#t)) #t))
(values Q (matrix-upper-triangle (matrix* (matrix-hermitian Q) M)))) (values Q (matrix-upper-triangle (matrix* (matrix-hermitian Q) M) zero)))
(: matrix-qr (case-> ((Matrix Real) -> (Values (Matrix Real) (Matrix Real))) (: matrix-qr (case-> ((Matrix Flonum) -> (Values (Matrix Flonum) (Matrix Flonum)))
((Matrix Flonum) Any -> (Values (Matrix Flonum) (Matrix Flonum)))
((Matrix Real) -> (Values (Matrix Real) (Matrix Real)))
((Matrix Real) Any -> (Values (Matrix Real) (Matrix Real))) ((Matrix Real) Any -> (Values (Matrix Real) (Matrix Real)))
((Matrix Float-Complex) -> (Values (Matrix Float-Complex)
(Matrix Float-Complex)))
((Matrix Float-Complex) Any -> (Values (Matrix Float-Complex)
(Matrix Float-Complex)))
((Matrix Number) -> (Values (Matrix Number) (Matrix Number))) ((Matrix Number) -> (Values (Matrix Number) (Matrix Number)))
((Matrix Number) Any -> (Values (Matrix Number) (Matrix Number))))) ((Matrix Number) Any -> (Values (Matrix Number) (Matrix Number)))))
(define (matrix-qr M [full? #t]) (define (matrix-qr M [full? #t])

38
collects/math/private/matrix/matrix-solve.rkt
View File

@ -24,7 +24,9 @@
;; =================================================================================================== ;; ===================================================================================================
;; Determinant ;; Determinant
(: matrix-determinant (case-> ((Matrix Real) -> Real) (: matrix-determinant (case-> ((Matrix Flonum) -> Flonum)
((Matrix Real) -> Real)
((Matrix Float-Complex) -> Float-Complex)
((Matrix Number) -> Number))) ((Matrix Number) -> Number)))
(define (matrix-determinant M) (define (matrix-determinant M)
(define m (square-matrix-size M)) (define m (square-matrix-size M))
@ -41,20 +43,22 @@
[else [else
(matrix-determinant/row-reduction M)])) (matrix-determinant/row-reduction M)]))
(: matrix-determinant/row-reduction (case-> ((Matrix Real) -> Real) (: matrix-determinant/row-reduction (case-> ((Matrix Flonum) -> Flonum)
((Matrix Real) -> Real)
((Matrix Float-Complex) -> Float-Complex)
((Matrix Number) -> Number))) ((Matrix Number) -> Number)))
(define (matrix-determinant/row-reduction M) (define (matrix-determinant/row-reduction M)
(define m (square-matrix-size M)) (define m (square-matrix-size M))
(define rows (matrix->vector* M)) (define rows (matrix->vector* M))
(let loop ([#{i : Nonnegative-Fixnum} 0] [#{sign : Real} 1]) (let loop ([#{i : Nonnegative-Fixnum} 0] [#{sign : (U Positive-Fixnum Negative-Fixnum)} 1])
(cond (cond
[(i . fx< . m) [(i . fx< . m)
(define-values (p pivot) (find-partial-pivot rows m i i)) (define-values (p pivot) (find-partial-pivot rows m i i))
(cond (cond
[(zero? pivot) 0] ; no pivot means non-invertible matrix [(zero? pivot) pivot] ; no pivot means non-invertible matrix
[else [else
(let ([sign (if (= i p) sign (begin (vector-swap! rows i p) ; swapping negates sign (let ([sign (if (= i p) sign (begin (vector-swap! rows i p) ; swapping negates sign
(* -1 sign)))]) (if (= sign 1) -1 1)))])
(elim-rows! rows m i i pivot (fx+ i 1)) ; adding scaled rows doesn't change it (elim-rows! rows m i i pivot (fx+ i 1)) ; adding scaled rows doesn't change it
(loop (fx+ i 1) sign))])] (loop (fx+ i 1) sign))])]
[else [else
@ -72,8 +76,12 @@
(and (square-matrix? M) (and (square-matrix? M)
(not (zero? (matrix-determinant M))))) (not (zero? (matrix-determinant M)))))
(: matrix-inverse (All (A) (case-> ((Matrix Real) -> (Matrix Real)) (: matrix-inverse (All (A) (case-> ((Matrix Flonum) -> (Matrix Flonum))
((Matrix Flonum) (-> A) -> (U A (Matrix Flonum)))
((Matrix Real) -> (Matrix Real))
((Matrix Real) (-> A) -> (U A (Matrix Real))) ((Matrix Real) (-> A) -> (U A (Matrix Real)))
((Matrix Float-Complex) -> (Matrix Float-Complex))
((Matrix Float-Complex) (-> A) -> (U A (Matrix Float-Complex)))
((Matrix Number) -> (Matrix Number)) ((Matrix Number) -> (Matrix Number))
((Matrix Number) (-> A) -> (U A (Matrix Number)))))) ((Matrix Number) (-> A) -> (U A (Matrix Number))))))
(define matrix-inverse (define matrix-inverse
@ -81,7 +89,8 @@
[(M) (matrix-inverse M (λ () (raise-argument-error 'matrix-inverse "matrix-invertible?" M)))] [(M) (matrix-inverse M (λ () (raise-argument-error 'matrix-inverse "matrix-invertible?" M)))]
[(M fail) [(M fail)
(define m (square-matrix-size M)) (define m (square-matrix-size M))
(define I (identity-matrix m)) (define x00 (matrix-ref M 0 0))
(define I (identity-matrix m (one* x00) (zero* x00)))
(define-values (IM^-1 wps) (parameterize ([array-strictness #f]) (define-values (IM^-1 wps) (parameterize ([array-strictness #f])
(matrix-gauss-elim (matrix-augment (list M I)) #t #t))) (matrix-gauss-elim (matrix-augment (list M I)) #t #t)))
(cond [(and (not (empty? wps)) (= (first wps) m)) (cond [(and (not (empty? wps)) (= (first wps) m))
@ -91,11 +100,16 @@
;; =================================================================================================== ;; ===================================================================================================
;; Solving linear systems ;; Solving linear systems
(: matrix-solve (All (A) (case-> (: matrix-solve
((Matrix Real) (Matrix Real) -> (Matrix Real)) (All (A) (case->
((Matrix Real) (Matrix Real) (-> A) -> (U A (Matrix Real))) ((Matrix Flonum) (Matrix Flonum) -> (Matrix Flonum))
((Matrix Number) (Matrix Number) -> (Matrix Number)) ((Matrix Flonum) (Matrix Flonum) (-> A) -> (U A (Matrix Flonum)))
((Matrix Number) (Matrix Number) (-> A) -> (U A (Matrix Number)))))) ((Matrix Real) (Matrix Real) -> (Matrix Real))
((Matrix Real) (Matrix Real) (-> A) -> (U A (Matrix Real)))
((Matrix Float-Complex) (Matrix Float-Complex) -> (Matrix Float-Complex))
((Matrix Float-Complex) (Matrix Float-Complex) (-> A) -> (U A (Matrix Float-Complex)))
((Matrix Number) (Matrix Number) -> (Matrix Number))
((Matrix Number) (Matrix Number) (-> A) -> (U A (Matrix Number))))))
(define matrix-solve (define matrix-solve
(case-lambda (case-lambda
[(M B) (matrix-solve M B (λ () (raise-argument-error 'matrix-solve "matrix-invertible?" 0 M B)))] [(M B) (matrix-solve M B (λ () (raise-argument-error 'matrix-solve "matrix-invertible?" 0 M B)))]

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

@ -35,7 +35,9 @@
(cond [(= j0 j1) Bs] (cond [(= j0 j1) Bs]
[else (cons (submatrix M (::) (:: j0 j1)) Bs)])) [else (cons (submatrix M (::) (:: j0 j1)) Bs)]))
(: matrix-col-space/ns (All (A) (case-> ((Matrix Real) -> (U #f (Matrix Real))) (: matrix-col-space/ns (All (A) (case-> ((Matrix Flonum) -> (U #f (Matrix Flonum)))
((Matrix Real) -> (U #f (Matrix Real)))
((Matrix Float-Complex) -> (U #f (Matrix Float-Complex)))
((Matrix Number) -> (U #f (Matrix Number)))))) ((Matrix Number) -> (U #f (Matrix Number))))))
(define (matrix-col-space/ns M) (define (matrix-col-space/ns M)
(define n (matrix-num-cols M)) (define n (matrix-num-cols M))
@ -52,13 +54,19 @@
(define next-j (first wps)) (define next-j (first wps))
(loop next-j (rest wps) (maybe-cons-submatrix M (fx+ j 1) next-j Bs))]))])) (loop next-j (rest wps) (maybe-cons-submatrix M (fx+ j 1) next-j Bs))]))]))
(: matrix-col-space (All (A) (case-> ((Matrix Real) -> (Matrix Real)) (: matrix-col-space (All (A) (case-> ((Matrix Flonum) -> (Array Flonum))
((Matrix Flonum) (-> A) -> (U A (Matrix Flonum)))
((Matrix Real) -> (Array Real))
((Matrix Real) (-> A) -> (U A (Matrix Real))) ((Matrix Real) (-> A) -> (U A (Matrix Real)))
((Matrix Number) -> (Matrix Number)) ((Matrix Float-Complex) -> (Array Float-Complex))
((Matrix Float-Complex) (-> A) -> (U A (Matrix Float-Complex)))
((Matrix Number) -> (Array Number))
((Matrix Number) (-> A) -> (U A (Matrix Number)))))) ((Matrix Number) (-> A) -> (U A (Matrix Number))))))
(define matrix-col-space (define matrix-col-space
(case-lambda (case-lambda
[(M) (matrix-col-space M (λ () (make-array (vector 0 (matrix-num-cols M)) 0)))] [(M) (matrix-col-space M (λ ()
(define x00 (matrix-ref M 0 0))
(make-array (vector 0 (matrix-num-cols M)) (zero* x00))))]
[(M fail) [(M fail)
(define S (parameterize ([array-strictness #f]) (define S (parameterize ([array-strictness #f])
(matrix-col-space/ns M))) (matrix-col-space/ns M)))

43
collects/math/private/matrix/typed-matrix-arithmetic.rkt
View File

@ -65,46 +65,65 @@
(loop (first arrs) (rest arrs)))])))])) (loop (first arrs) (rest arrs)))])))]))
(: matrix*/ns (case-> ((Matrix Real) (Listof (Matrix Real)) -> (Matrix Real)) (: matrix*/ns
((Matrix Number) (Listof (Matrix Number)) -> (Matrix Number)))) (case-> ((Matrix Flonum) (Listof (Matrix Flonum)) -> (Matrix Flonum))
((Matrix Real) (Listof (Matrix Real)) -> (Matrix Real))
((Matrix Float-Complex) (Listof (Matrix Float-Complex)) -> (Matrix Float-Complex))
((Matrix Number) (Listof (Matrix Number)) -> (Matrix Number))))
(define (matrix*/ns a as) (define (matrix*/ns a as)
(cond [(empty? as) a] (cond [(empty? as) a]
[else (matrix*/ns (inline-matrix-multiply a (first as)) (rest as))])) [else (matrix*/ns (inline-matrix-multiply a (first as)) (rest as))]))
(: matrix* (case-> ((Matrix Real) (Matrix Real) * -> (Matrix Real)) (: matrix* (case-> ((Matrix Flonum) (Matrix Flonum) * -> (Matrix Flonum))
((Matrix Real) (Matrix Real) * -> (Matrix Real))
((Matrix Float-Complex) (Matrix Float-Complex) * -> (Matrix Float-Complex))
((Matrix Number) (Matrix Number) * -> (Matrix Number)))) ((Matrix Number) (Matrix Number) * -> (Matrix Number))))
(define (matrix* a . as) (call/ns (λ () (matrix*/ns a as)))) (define (matrix* a . as) (call/ns (λ () (matrix*/ns a as))))
(: matrix+/ns (case-> ((Matrix Real) (Listof (Matrix Real)) -> (Matrix Real)) (: matrix+/ns
((Matrix Number) (Listof (Matrix Number)) -> (Matrix Number)))) (case-> ((Matrix Flonum) (Listof (Matrix Flonum)) -> (Matrix Flonum))
((Matrix Real) (Listof (Matrix Real)) -> (Matrix Real))
((Matrix Float-Complex) (Listof (Matrix Float-Complex)) -> (Matrix Float-Complex))
((Matrix Number) (Listof (Matrix Number)) -> (Matrix Number))))
(define (matrix+/ns a as) (define (matrix+/ns a as)
(cond [(empty? as) a] (cond [(empty? as) a]
[else (matrix+/ns (inline-matrix+ a (first as)) (rest as))])) [else (matrix+/ns (inline-matrix+ a (first as)) (rest as))]))
(: matrix+ (case-> ((Matrix Real) (Matrix Real) * -> (Matrix Real)) (: matrix+ (case-> ((Matrix Flonum) (Matrix Flonum) * -> (Matrix Flonum))
((Matrix Real) (Matrix Real) * -> (Matrix Real))
((Matrix Float-Complex) (Matrix Float-Complex) * -> (Matrix Float-Complex))
((Matrix Number) (Matrix Number) * -> (Matrix Number)))) ((Matrix Number) (Matrix Number) * -> (Matrix Number))))
(define (matrix+ a . as) (call/ns (λ () (matrix+/ns a as)))) (define (matrix+ a . as) (call/ns (λ () (matrix+/ns a as))))
(: matrix-/ns (case-> ((Matrix Real) (Listof (Matrix Real)) -> (Matrix Real)) (: matrix-/ns
((Matrix Number) (Listof (Matrix Number)) -> (Matrix Number)))) (case-> ((Matrix Flonum) (Listof (Matrix Flonum)) -> (Matrix Flonum))
((Matrix Real) (Listof (Matrix Real)) -> (Matrix Real))
((Matrix Float-Complex) (Listof (Matrix Float-Complex)) -> (Matrix Float-Complex))
((Matrix Number) (Listof (Matrix Number)) -> (Matrix Number))))
(define (matrix-/ns a as) (define (matrix-/ns a as)
(cond [(empty? as) a] (cond [(empty? as) a]
[else (matrix-/ns (inline-matrix- a (first as)) (rest as))])) [else (matrix-/ns (inline-matrix- a (first as)) (rest as))]))
(: matrix- (case-> ((Matrix Real) (Matrix Real) * -> (Matrix Real)) (: matrix- (case-> ((Matrix Flonum) (Matrix Flonum) * -> (Matrix Flonum))
((Matrix Real) (Matrix Real) * -> (Matrix Real))
((Matrix Float-Complex) (Matrix Float-Complex) * -> (Matrix Float-Complex))
((Matrix Number) (Matrix Number) * -> (Matrix Number)))) ((Matrix Number) (Matrix Number) * -> (Matrix Number))))
(define (matrix- a . as) (define (matrix- a . as)
(call/ns (λ () (cond [(empty? as) (inline-matrix- a)] (call/ns (λ () (cond [(empty? as) (inline-matrix-scale a -1)]
[else (matrix-/ns a as)])))) [else (matrix-/ns a as)]))))
(: matrix-scale (case-> ((Matrix Real) Real -> (Matrix Real)) (: matrix-scale (case-> ((Matrix Flonum) Flonum -> (Matrix Flonum))
((Matrix Real) Real -> (Matrix Real))
((Matrix Float-Complex) Float-Complex -> (Matrix Float-Complex))
((Matrix Number) Number -> (Matrix Number)))) ((Matrix Number) Number -> (Matrix Number))))
(define (matrix-scale a x) (inline-matrix-scale a x)) (define (matrix-scale a x) (inline-matrix-scale a x))
(: matrix-sum (case-> ((Listof (Matrix Real)) -> (Matrix Real)) (: matrix-sum (case-> ((Listof (Matrix Flonum)) -> (Matrix Flonum))
((Listof (Matrix Real)) -> (Matrix Real))
((Listof (Matrix Float-Complex)) -> (Matrix Float-Complex))
((Listof (Matrix Number)) -> (Matrix Number)))) ((Listof (Matrix Number)) -> (Matrix Number))))
(define (matrix-sum lst) (define (matrix-sum lst)
(cond [(empty? lst) (raise-argument-error 'matrix-sum "nonempty List" lst)] (cond [(empty? lst) (raise-argument-error 'matrix-sum "nonempty List" lst)]

92
collects/math/private/matrix/untyped-matrix-arithmetic.rkt
View File

@ -8,39 +8,67 @@
inline-matrix-map inline-matrix-map
matrix-map) matrix-map)
(module syntax-defs racket/base (module typed-multiply-defs typed/racket/base
(require (for-syntax racket/base) (require racket/fixnum
(only-in typed/racket/base λ: : inst Index)
"matrix-types.rkt" "matrix-types.rkt"
"utils.rkt" "utils.rkt"
"../unsafe.rkt"
"../array/array-struct.rkt" "../array/array-struct.rkt"
"../array/array-fold.rkt"
"../array/array-transform.rkt" "../array/array-transform.rkt"
"../array/mutable-array.rkt"
"../array/utils.rkt") "../array/utils.rkt")
(provide (all-defined-out)) (provide (all-defined-out))
;(: matrix-multiply ((Matrix Number) (Matrix Number) -> (Matrix Number))) (: matrix-multiply-data (All (A) ((Matrix A) (Matrix A) -> (Values Index Index Index
(Vectorof A) (Vectorof A)
(-> (Boxof A))))))
(define (matrix-multiply-data arr brr)
(let-values ([(m p n) (matrix-multiply-shape arr brr)])
(define arr-data (mutable-array-data (array->mutable-array arr)))
(define brr-data (mutable-array-data (array->mutable-array (parameterize ([array-strictness #f])
(array-axis-swap brr 0 1)))))
(define bx (make-thread-local-box (unsafe-vector-ref arr-data 0)))
(values m p n arr-data brr-data bx)))
(: make-matrix-multiply (All (A) (Index Index Index (Index Index -> A) -> (Matrix A))))
(define (make-matrix-multiply m p n sum-loop)
(array-default-strict
(unsafe-build-array
((inst vector Index) m n)
(λ: ([ij : Indexes])
(sum-loop (assert (fx* (unsafe-vector-ref ij 0) p) index?)
(assert (fx* (unsafe-vector-ref ij 1) p) index?))))))
) ; module
(module untyped-multiply-defs racket/base
(require (for-syntax racket/base)
racket/fixnum
racket/unsafe/ops
(only-in typed/racket/base λ: : Index let: Nonnegative-Fixnum)
(submod ".." typed-multiply-defs)
"matrix-types.rkt"
"utils.rkt"
"../array/array-struct.rkt")
(provide (all-defined-out))
;; This is a macro so the result can have as precise a type as possible ;; This is a macro so the result can have as precise a type as possible
(define-syntax-rule (inline-matrix-multiply arr-expr brr-expr) (define-syntax-rule (inline-matrix-multiply arr brr)
(let ([arr arr-expr] (let-values ([(m p n arr-data brr-data bx) (matrix-multiply-data arr brr)])
[brr brr-expr]) (make-matrix-multiply
(let-values ([(m p n) (matrix-multiply-shape arr brr)] m p n
;; Make arr strict because its elements are reffed multiple times (λ: ([i : Index] [j : Index])
[(_) (array-strict! arr)]) (let ([bx (bx)]
(let (;; Extend arr in the center dimension [v (* (unsafe-vector-ref arr-data i)
[arr-proc (unsafe-array-proc (array-axis-insert arr 1 n))] (unsafe-vector-ref brr-data j))])
;; Transpose brr and extend in the leftmost dimension (let: loop ([k : Nonnegative-Fixnum 1] [v v])
[brr-proc (unsafe-array-proc (cond [(k . fx< . p)
(array-axis-insert (array-strict (array-axis-swap brr 0 1)) 0 m))]) (loop (fx+ k 1)
;; The *transpose* of brr is traversed in row-major order when this result is traversed (+ v (* (unsafe-vector-ref arr-data (fx+ i k))
;; in row-major order (which is why the transpose is strict, not brr) (unsafe-vector-ref brr-data (fx+ j k)))))]
(array-axis-sum [else (set-box! bx v)]))
(unsafe-build-array (unbox bx))))))
((inst vector Index) m n p)
(λ: ([js : Indexes])
(* (arr-proc js) (brr-proc js))))
2)))))
(define-syntax (do-inline-matrix* stx) (define-syntax (do-inline-matrix* stx)
(syntax-case stx () (syntax-case stx ()
@ -54,6 +82,19 @@
(parameterize ([array-strictness #f]) (parameterize ([array-strictness #f])
(do-inline-matrix* arr brrs ...)))) (do-inline-matrix* arr brrs ...))))
) ; module
(module syntax-defs racket/base
(require (for-syntax racket/base)
(only-in typed/racket/base λ: : inst Index)
(submod ".." typed-multiply-defs)
"matrix-types.rkt"
"utils.rkt"
"../array/array-struct.rkt"
"../array/utils.rkt")
(provide (all-defined-out))
(define-syntax (inline-matrix-map stx) (define-syntax (inline-matrix-map stx)
(syntax-case stx () (syntax-case stx ()
[(_ f arr-expr) [(_ f arr-expr)
@ -117,5 +158,6 @@
) ; module ) ; module
(require 'syntax-defs (require 'untyped-multiply-defs
'syntax-defs
'untyped-defs) 'untyped-defs)

66
collects/math/private/matrix/utils.rkt
View File

@ -3,6 +3,7 @@
(require racket/performance-hint (require racket/performance-hint
racket/string racket/string
racket/fixnum racket/fixnum
math/base
"matrix-types.rkt" "matrix-types.rkt"
"../unsafe.rkt" "../unsafe.rkt"
"../array/array-struct.rkt" "../array/array-struct.rkt"
@ -69,7 +70,9 @@
(rational? (imag-part z)))]))) (rational? (imag-part z)))])))
(: find-partial-pivot (: find-partial-pivot
(case-> ((Vectorof (Vectorof Real)) Index Index Index -> (Values Index Real)) (case-> ((Vectorof (Vectorof Flonum)) Index Index Index -> (Values Index Flonum))
((Vectorof (Vectorof Real)) Index Index Index -> (Values Index Real))
((Vectorof (Vectorof Float-Complex)) Index Index Index -> (Values Index Float-Complex))
((Vectorof (Vectorof Number)) Index Index Index -> (Values Index Number)))) ((Vectorof (Vectorof Number)) Index Index Index -> (Values Index Number))))
;; Find the element with maximum magnitude in a column ;; Find the element with maximum magnitude in a column
(define (find-partial-pivot rows m i j) (define (find-partial-pivot rows m i j)
@ -85,7 +88,9 @@
[else (values p pivot)]))) [else (values p pivot)])))
(: find-first-pivot (: find-first-pivot
(case-> ((Vectorof (Vectorof Real)) Index Index Index -> (Values Index Real)) (case-> ((Vectorof (Vectorof Flonum)) Index Index Index -> (Values Index Flonum))
((Vectorof (Vectorof Real)) Index Index Index -> (Values Index Real))
((Vectorof (Vectorof Float-Complex)) Index Index Index -> (Values Index Float-Complex))
((Vectorof (Vectorof Number)) Index Index Index -> (Values Index Number)))) ((Vectorof (Vectorof Number)) Index Index Index -> (Values Index Number))))
;; Find the first nonzero element in a column ;; Find the first nonzero element in a column
(define (find-first-pivot rows m i j) (define (find-first-pivot rows m i j)
@ -100,7 +105,10 @@
(values i pivot)])))) (values i pivot)]))))
(: elim-rows! (: elim-rows!
(case-> ((Vectorof (Vectorof Real)) Index Index Index Real Nonnegative-Fixnum -> Void) (case-> ((Vectorof (Vectorof Flonum)) Index Index Index Flonum Nonnegative-Fixnum -> Void)
((Vectorof (Vectorof Real)) Index Index Index Real Nonnegative-Fixnum -> Void)
((Vectorof (Vectorof Float-Complex)) Index Index Index Float-Complex Nonnegative-Fixnum
-> Void)
((Vectorof (Vectorof Number)) Index Index Index Number Nonnegative-Fixnum -> Void))) ((Vectorof (Vectorof Number)) Index Index Index Number Nonnegative-Fixnum -> Void)))
(define (elim-rows! rows m i j pivot start) (define (elim-rows! rows m i j pivot start)
(define row_i (unsafe-vector-ref rows i)) (define row_i (unsafe-vector-ref rows i))
@ -109,8 +117,8 @@
(unless (l . fx= . i) (unless (l . fx= . i)
(define row_l (unsafe-vector-ref rows l)) (define row_l (unsafe-vector-ref rows l))
(define x_lj (unsafe-vector-ref row_l j)) (define x_lj (unsafe-vector-ref row_l j))
(unless (zero? x_lj) (unless (= x_lj 0)
(vector-scaled-add! row_l row_i (- (/ x_lj pivot)) j) (vector-scaled-add! row_l row_i (* -1 (/ x_lj pivot)) j)
;; Make sure the element below the pivot is zero ;; Make sure the element below the pivot is zero
(unsafe-vector-set! row_l j (- x_lj x_lj)))) (unsafe-vector-set! row_l j (- x_lj x_lj))))
(loop (fx+ l 1))))) (loop (fx+ l 1)))))
@ -124,3 +132,51 @@
(thnk)))) (thnk))))
) ; begin-encourage-inline ) ; begin-encourage-inline
(: make-thread-local-box (All (A) (A -> (-> (Boxof A)))))
(define (make-thread-local-box contents)
(let: ([val : (Thread-Cellof (U #f (Boxof A))) (make-thread-cell #f)])
(λ () (or (thread-cell-ref val)
(let: ([v : (Boxof A) (box contents)])
(thread-cell-set! val v)
v)))))
(: one (case-> (Flonum -> Nonnegative-Flonum)
(Real -> (U 1 Nonnegative-Flonum))
(Float-Complex -> Nonnegative-Flonum)
(Number -> (U 1 Nonnegative-Flonum))))
(define (one x)
(cond [(flonum? x) 1.0]
[(real? x) 1]
[(float-complex? x) 1.0]
[else 1]))
(: zero (case-> (Flonum -> Flonum-Positive-Zero)
(Real -> (U 0 Flonum-Positive-Zero))
(Float-Complex -> Flonum-Positive-Zero)
(Number -> (U 0 Flonum-Positive-Zero))))
(define (zero x)
(cond [(flonum? x) 0.0]
[(real? x) 0]
[(float-complex? x) 0.0]
[else 0]))
(: one* (case-> (Flonum -> Nonnegative-Flonum)
(Real -> (U 1 Nonnegative-Flonum))
(Float-Complex -> Float-Complex)
(Number -> (U 1 Nonnegative-Flonum Float-Complex))))
(define (one* x)
(cond [(flonum? x) 1.0]
[(real? x) 1]
[(float-complex? x) 1.0+0.0i]
[else 1]))
(: zero* (case-> (Flonum -> Flonum-Positive-Zero)
(Real -> (U 0 Flonum-Positive-Zero))
(Float-Complex -> Float-Complex)
(Number -> (U 0 Flonum-Positive-Zero Float-Complex))))
(define (zero* x)
(cond [(flonum? x) 0.0]
[(real? x) 0]
[(float-complex? x) 0.0+0.0i]
[else 0]))

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

@ -1,7 +1,7 @@
#lang typed/racket/base #lang typed/racket/base
(require racket/fixnum (require racket/fixnum
racket/math math/base
math/private/unsafe) math/private/unsafe)
(provide vector-swap! (provide vector-swap!
@ -14,9 +14,12 @@
vector-zero! vector-zero!
vector-zero?) vector-zero?)
(: mag^2 (Number -> Nonnegative-Real)) (: mag^2 (case-> (Flonum -> Nonnegative-Flonum)
(Float-Complex -> Nonnegative-Flonum)
(Number -> Nonnegative-Real)))
(define (mag^2 x) (define (mag^2 x)
(max 0 (real-part (* x (conjugate x))))) (cond [(real? x) (sqr x)]
[else (abs (real-part (* x (conjugate x))))]))
(: vector-swap! (All (A) ((Vectorof A) Integer Integer -> Void))) (: vector-swap! (All (A) ((Vectorof A) Integer Integer -> Void)))
(define (vector-swap! vs i0 i1) (define (vector-swap! vs i0 i1)
@ -35,7 +38,9 @@
(loop (fx+ i 1))) (loop (fx+ i 1)))
(void))))) (void)))))
(: vector-scale! (case-> ((Vectorof Real) Real -> Void) (: vector-scale! (case-> ((Vectorof Flonum) Flonum -> Void)
((Vectorof Real) Real -> Void)
((Vectorof Float-Complex) Float-Complex -> Void)
((Vectorof Number) Number -> Void))) ((Vectorof Number) Number -> Void)))
(define (vector-scale! vs v) (define (vector-scale! vs v)
(vector-generic-scale! vs v *)) (vector-generic-scale! vs v *))
@ -52,25 +57,38 @@
(loop (fx+ i 1))) (loop (fx+ i 1)))
(void))))) (void)))))
(: vector-scaled-add! (case-> ((Vectorof Real) (Vectorof Real) Real -> Void) (: vector-scaled-add!
((Vectorof Real) (Vectorof Real) Real Index -> Void) (case-> ((Vectorof Flonum) (Vectorof Flonum) Flonum -> Void)
((Vectorof Number) (Vectorof Number) Number -> Void) ((Vectorof Flonum) (Vectorof Flonum) Flonum Index -> Void)
((Vectorof Number) (Vectorof Number) Number Index -> Void))) ((Vectorof Real) (Vectorof Real) Real -> Void)
((Vectorof Real) (Vectorof Real) Real Index -> Void)
((Vectorof Float-Complex) (Vectorof Float-Complex) Float-Complex -> Void)
((Vectorof Float-Complex) (Vectorof Float-Complex) Float-Complex Index -> Void)
((Vectorof Number) (Vectorof Number) Number -> Void)
((Vectorof Number) (Vectorof Number) Number Index -> Void)))
(define (vector-scaled-add! vs0 vs1 s [start 0]) (define (vector-scaled-add! vs0 vs1 s [start 0])
(vector-generic-scaled-add! vs0 vs1 s start + *)) (vector-generic-scaled-add! vs0 vs1 s start + *))
(: vector-mag^2 (case-> ((Vectorof Real) -> Nonnegative-Real) (: vector-mag^2 (case-> ((Vectorof Flonum) -> Nonnegative-Flonum)
((Vectorof Real) -> Nonnegative-Real)
((Vectorof Float-Complex) -> Nonnegative-Flonum)
((Vectorof Number) -> Nonnegative-Real))) ((Vectorof Number) -> Nonnegative-Real)))
(define (vector-mag^2 vs) (define (vector-mag^2 vs)
(define n (vector-length vs)) (define n (vector-length vs))
(let loop ([#{i : Nonnegative-Fixnum} 0] [#{s : Nonnegative-Real} 0]) (cond [(fx= n 0) (raise-argument-error 'vector-mag^2 "nonempty Vector" vs)]
(if (i . fx>= . n) s (loop (fx+ i 1) (+ s (mag^2 (unsafe-vector-ref vs i))))))) [else
(define s (mag^2 (unsafe-vector-ref vs 0)))
(let: loop ([i : Nonnegative-Fixnum 1] [s s])
(cond [(i . fx< . n) (loop (fx+ i 1) (+ s (mag^2 (unsafe-vector-ref vs i))))]
[else (abs s)]))]))
(: vector-dot (case-> ((Vectorof Real) (Vectorof Real) -> Real) (: vector-dot (case-> ((Vectorof Flonum) (Vectorof Flonum) -> Flonum)
((Vectorof Real) (Vectorof Real) -> Real)
((Vectorof Float-Complex) (Vectorof Float-Complex) -> Float-Complex)
((Vectorof Number) (Vectorof Number) -> Number))) ((Vectorof Number) (Vectorof Number) -> Number)))
(define (vector-dot vs0 vs1) (define (vector-dot vs0 vs1)
(define n (min (vector-length vs0) (vector-length vs1))) (define n (min (vector-length vs0) (vector-length vs1)))
(cond [(= n 0) 0] (cond [(fx= n 0) (raise-argument-error 'vector-dot "nonempty Vector" 0 vs0 vs1)]
[else [else
(define v0 (unsafe-vector-ref vs0 0)) (define v0 (unsafe-vector-ref vs0 0))
(define v1 (unsafe-vector-ref vs1 0)) (define v1 (unsafe-vector-ref vs1 0))
@ -81,7 +99,9 @@
(loop (fx+ i 1) (+ s (* v0 (conjugate v1))))] (loop (fx+ i 1) (+ s (* v0 (conjugate v1))))]
[else s]))])) [else s]))]))
(: vector-normalize! (case-> ((Vectorof Real) -> Nonnegative-Real) (: vector-normalize! (case-> ((Vectorof Flonum) -> Nonnegative-Flonum)
((Vectorof Real) -> Nonnegative-Real)
((Vectorof Float-Complex) -> Nonnegative-Flonum)
((Vectorof Number) -> Nonnegative-Real))) ((Vectorof Number) -> Nonnegative-Real)))
(define (vector-normalize! vs) (define (vector-normalize! vs)
(define n (vector-length vs)) (define n (vector-length vs))
@ -93,31 +113,51 @@
(loop (fx+ i 1))))) (loop (fx+ i 1)))))
s) s)
(: vector-sub-proj! (case-> ((Vectorof Real) (Vectorof Real) Any -> Nonnegative-Real) (: one (case-> (Flonum -> Nonnegative-Flonum)
((Vectorof Number) (Vectorof Number) Any -> Nonnegative-Real))) (Real -> Nonnegative-Real)
(Float-Complex -> Nonnegative-Flonum)
(Number -> Nonnegative-Real)))
(define (one x)
(cond [(flonum? x) 1.0]
[(real? x) 1]
[(float-complex? x) 1.0]
[else 1]))
(: vector-sub-proj!
(case-> ((Vectorof Flonum) (Vectorof Flonum) Any -> Nonnegative-Flonum)
((Vectorof Real) (Vectorof Real) Any -> Nonnegative-Real)
((Vectorof Float-Complex) (Vectorof Float-Complex) Any -> Nonnegative-Flonum)
((Vectorof Number) (Vectorof Number) Any -> Nonnegative-Real)))
(define (vector-sub-proj! vs0 vs1 unit?) (define (vector-sub-proj! vs0 vs1 unit?)
(define n (min (vector-length vs0) (vector-length vs1))) (define n (min (vector-length vs0) (vector-length vs1)))
(define t (if unit? 1 (vector-mag^2 vs1))) (cond [(fx= n 0) (raise-argument-error 'vector-sub-proj! "nonempty Vector" 0 vs0 vs1)]
(unless (and (zero? t) (exact? t)) [else
(define s (/ (vector-dot vs0 vs1) t)) (define t (if unit? (one (unsafe-vector-ref vs0 0)) (vector-mag^2 vs1)))
(let loop ([#{i : Nonnegative-Fixnum} 0]) (unless (and (zero? t) (exact? t))
(when (i . fx< . n) (define s (/ (vector-dot vs0 vs1) t))
(define v0 (unsafe-vector-ref vs0 i)) (let loop ([#{i : Nonnegative-Fixnum} 0])
(define v1 (unsafe-vector-ref vs1 i)) (when (i . fx< . n)
(unsafe-vector-set! vs0 i (- v0 (* v1 s))) (define v0 (unsafe-vector-ref vs0 i))
(loop (fx+ i 1))))) (define v1 (unsafe-vector-ref vs1 i))
t) (unsafe-vector-set! vs0 i (- v0 (* v1 s)))
(loop (fx+ i 1)))))
t]))
(: vector-zero! (case-> ((Vectorof Real) -> Void) (: vector-zero! (case-> ((Vectorof Flonum) -> Void)
((Vectorof Real) -> Void)
((Vectorof Float-Complex) -> Void)
((Vectorof Number) -> Void))) ((Vectorof Number) -> Void)))
(define (vector-zero! vs) (define (vector-zero! vs)
(define n (vector-length vs)) (define n (vector-length vs))
(let loop ([#{i : Nonnegative-Fixnum} 0]) (let loop ([#{i : Nonnegative-Fixnum} 0])
(when (i . fx< . n) (when (i . fx< . n)
(unsafe-vector-set! vs i 0) (define x (unsafe-vector-ref vs i))
(unsafe-vector-set! vs i (- x x))
(loop (fx+ i 1))))) (loop (fx+ i 1)))))
(: vector-zero? (case-> ((Vectorof Real) -> Boolean) (: vector-zero? (case-> ((Vectorof Flonum) -> Boolean)
((Vectorof Real) -> Boolean)
((Vectorof Float-Complex) -> Boolean)
((Vectorof Number) -> Boolean))) ((Vectorof Number) -> Boolean)))
(define (vector-zero? vs) (define (vector-zero? vs)
(define n (vector-length vs)) (define n (vector-length vs))

76
collects/math/scribblings/math-matrix.scrbl
View File

@ -45,12 +45,13 @@ Like all of @racketmodname[math], @racketmodname[math/matrix] is a work in progr
Most of the basic algorithms are implemented, but some are still in planning. Most of the basic algorithms are implemented, but some are still in planning.
Possibly the most useful unimplemented algorithms are Possibly the most useful unimplemented algorithms are
@itemlist[@item{LUP decomposition (currently, LU decomposition is implemented, in @racket[matrix-lu])} @itemlist[@item{LUP decomposition (currently, LU decomposition is implemented, in @racket[matrix-lu])}
@item{@racket[matrix-solve] for upper-triangular matrices} @item{@racket[matrix-solve] for triangular matrices}
@item{Singular value decomposition (SVD)} @item{Singular value decomposition (SVD)}
@item{Eigendecomposition} @item{Eigendecomposition}
@item{Decomposition-based solvers} @item{Decomposition-based solvers}
@item{Pseudoinverse, least-squares fitting}] @item{Pseudoinverse and least-squares solving}]
All of these are planned for the next Racket release. All of these are planned for the next Racket release, as well as fast flonum-specific matrix
operations and LAPACK integration.
@local-table-of-contents[] @local-table-of-contents[]
@ -66,8 +67,7 @@ From the point of view of the functions in @racketmodname[math/matrix], a @defte
Technically, a matrix's entries may be any type, and some fully polymorphic matrix functions such as Technically, a matrix's entries may be any type, and some fully polymorphic matrix functions such as
@racket[matrix-row] and @racket[matrix-map] operate on any kind of matrix. @racket[matrix-row] and @racket[matrix-map] operate on any kind of matrix.
Other functions, such as @racket[matrix+], require their input matrices to contain either Other functions, such as @racket[matrix+], require their matrix arguments to contain numeric values.
@racket[Real] or @racket[Number] values.
@subsection[#:tag "matrix:function-types"]{Function Types} @subsection[#:tag "matrix:function-types"]{Function Types}
@ -77,17 +77,22 @@ that a return value is a matrix.
Most functions that implement matrix algorithms are documented as accepting @racket[(Matrix Number)] Most functions that implement matrix algorithms are documented as accepting @racket[(Matrix Number)]
values. This includes @racket[(Matrix Real)], which is a subtype. Most of these functions have a more values. This includes @racket[(Matrix Real)], which is a subtype. Most of these functions have a more
precise type than is documented. For example, @racket[matrix-conjugate] actually has the type precise type than is documented. For example, @racket[matrix-conjugate] has the type
@racketblock[(case-> ((Matrix Real) -> (Matrix Real)) @racketblock[(case-> ((Matrix Flonum) -> (Matrix Flonum))
((Matrix Number) -> (Matrix Number)))] ((Matrix Real) -> (Matrix Real))
even though it is documented as having the less precise type ((Matrix Float-Complex) -> (Matrix Float-Complex))
@racket[((Matrix Number) -> (Matrix Number))]. ((Matrix Number) -> (Matrix Number)))]
but is documented as having the type @racket[((Matrix Number) -> (Matrix Number))].
Precise function types allow Typed Racket to prove more facts about @racketmodname[math/matrix] Precise function types allow Typed Racket to prove more facts about @racketmodname[math/matrix]
client programs. In particular, it is usually easy for it to prove that matrix expressions on real client programs. In particular, it is usually easy for it to prove that operations on real
matrices return real matrices: matrices return real matrices:
@interaction[#:eval typed-eval @interaction[#:eval typed-eval
(matrix-conjugate (matrix [[1 2 3] [4 5 6]]))] (matrix-conjugate (matrix [[1 2 3] [4 5 6]]))]
and that operations on inexact matrices return inexact matrices:
@interaction[#:eval typed-eval
(matrix-conjugate (matrix [[1.0+2.0i 2.0+3.0i 3.0+4.0i]
[4.0+5.0i 5.0+6.0i 6.0+7.0i]]))]
@subsection[#:tag "matrix:failure"]{Failure Arguments} @subsection[#:tag "matrix:failure"]{Failure Arguments}
@ -102,7 +107,7 @@ For example, the (simplified) type of @racket[matrix-inverse] is
Thus, if a failure thunk is given, the call site is required to check for return values of type Thus, if a failure thunk is given, the call site is required to check for return values of type
@racket[F] explicitly. @racket[F] explicitly.
The default failure thunk, which raises an error, has type @racket[(-> Nothing)]. Default failure thunks usually raise an error, and have the type @racket[(-> Nothing)].
For such failure thunks, @racket[(U F (Matrix Number))] is equivalent to @racket[(Matrix Number)], For such failure thunks, @racket[(U F (Matrix Number))] is equivalent to @racket[(Matrix Number)],
because @racket[Nothing] is part of every type. (In Racket, any expression may raise an error.) because @racket[Nothing] is part of every type. (In Racket, any expression may raise an error.)
Thus, in this case, no explicit test for values of type @racket[F] is necessary (though of course they Thus, in this case, no explicit test for values of type @racket[F] is necessary (though of course they
@ -110,8 +115,8 @@ may be caught using @racket[with-handlers] or similar).
@subsection[#:tag "matrix:broadcasting"]{Broadcasting} @subsection[#:tag "matrix:broadcasting"]{Broadcasting}
Pointwise matrix operations do not @tech{broadcast} their arguments when given matrices with Unlike array operations, pointwise matrix operations @bold{do not} @tech{broadcast} their arguments
different sizes: when given matrices with different axis lengths:
@interaction[#:eval typed-eval @interaction[#:eval typed-eval
(matrix+ (identity-matrix 2) (matrix [[10]]))] (matrix+ (identity-matrix 2) (matrix [[10]]))]
If you need broadcasting, use array operations: If you need broadcasting, use array operations:
@ -217,8 +222,12 @@ Like @racket[matrix], but returns a @tech{column matrix}.
(col-matrix [])] (col-matrix [])]
} }
@defproc[(identity-matrix [n Integer]) (Matrix (U 0 1))]{ @defproc[(identity-matrix [n Integer] [one A 1] [zero A 0]) (Matrix A)]{
Returns an @racket[n]×@racket[n] identity matrix; @racket[n] must be positive. Returns an @racket[n]×@racket[n] identity matrix, which has the value @racket[one] on the diagonal
and @racket[zero] everywhere else. The height/width @racket[n] must be positive.
@examples[#:eval typed-eval
(identity-matrix 3)
(identity-matrix 4 1.0+0.0i 0.0+0.0i)]
} }
@defproc[(make-matrix [m Integer] [n Integer] [x A]) (Matrix A)]{ @defproc[(make-matrix [m Integer] [n Integer] [x A]) (Matrix A)]{
@ -233,22 +242,28 @@ both @racket[m] and @racket[n] must be positive.
Analogous to @racket[build-array] (and defined in terms of it). Analogous to @racket[build-array] (and defined in terms of it).
} }
@defproc[(diagonal-matrix [xs (Listof A)]) (Matrix (U A 0))]{ @defproc[(diagonal-matrix [xs (Listof A)] [zero A 0]) (Matrix A)]{
Returns a matrix with @racket[xs] along the diagonal and @racket[0] everywhere else. Returns a matrix with @racket[xs] along the diagonal and @racket[zero] everywhere else.
The length of @racket[xs] must be positive. The length of @racket[xs] must be positive.
@examples[#:eval typed-eval
(diagonal-matrix '(1 2 3 4 5 6))
(diagonal-matrix '(1.0 2.0 3.0 4.0 5.0) 0.0)]
} }
@define[block-diagonal-url]{http://en.wikipedia.org/wiki/Block_matrix#Block_diagonal_matrices} @define[block-diagonal-url]{http://en.wikipedia.org/wiki/Block_matrix#Block_diagonal_matrices}
@margin-note*{@hyperlink[block-diagonal-url]{Wikipedia: Block-diagonal matrices}} @margin-note*{@hyperlink[block-diagonal-url]{Wikipedia: Block-diagonal matrices}}
@defproc[(block-diagonal-matrix [Xs (Listof (Matrix A))]) (Matrix (U A 0))]{ @defproc[(block-diagonal-matrix [Xs (Listof (Matrix A))] [zero A 0]) (Matrix A)]{
Returns a matrix with matrices @racket[Xs] along the diagonal and @racket[0] everywhere else. Returns a matrix with matrices @racket[Xs] along the diagonal and @racket[zero] everywhere else.
The length of @racket[Xs] must be positive. The length of @racket[Xs] must be positive.
@examples[#:eval typed-eval @examples[#:eval typed-eval
(block-diagonal-matrix (list (matrix [[6 7] [8 9]]) (block-diagonal-matrix (list (matrix [[6 7] [8 9]])
(diagonal-matrix '(7 5 7)) (diagonal-matrix '(7 5 7))
(col-matrix [1 2 3]) (col-matrix [1 2 3])
(row-matrix [4 5 6])))] (row-matrix [4 5 6])))
(block-diagonal-matrix (list (make-matrix 2 2 2.0+3.0i)
(make-matrix 2 2 5.0+7.0i))
0.0+0.0i)]
} }
@define[vandermonde-url]{http://en.wikipedia.org/wiki/Vandermonde_matrix} @define[vandermonde-url]{http://en.wikipedia.org/wiki/Vandermonde_matrix}
@ -257,7 +272,8 @@ The length of @racket[Xs] must be positive.
@defproc[(vandermonde-matrix [xs (Listof Number)] [n Integer]) (Matrix Number)]{ @defproc[(vandermonde-matrix [xs (Listof Number)] [n Integer]) (Matrix Number)]{
Returns an @racket[m]×@racket[n] Vandermonde matrix, where @racket[m = (length xs)]. Returns an @racket[m]×@racket[n] Vandermonde matrix, where @racket[m = (length xs)].
@examples[#:eval typed-eval @examples[#:eval typed-eval
(vandermonde-matrix '(1 2 3 4) 5)] (vandermonde-matrix '(1 2 3 4) 5)
(vandermonde-matrix '(5.2 3.4 2.0) 3)]
Using a Vandermonde matrix to find a Lagrange polynomial (the polynomial of least degree that passes Using a Vandermonde matrix to find a Lagrange polynomial (the polynomial of least degree that passes
through a given set of points): through a given set of points):
@interaction[#:eval untyped-eval @interaction[#:eval untyped-eval
@ -384,7 +400,8 @@ Computes @racket[(matrix* M ...)] with @racket[n] arguments, but more efficientl
@examples[#:eval untyped-eval @examples[#:eval untyped-eval
; The 100th (and 101th) Fibonacci number: ; The 100th (and 101th) Fibonacci number:
(matrix* (matrix-expt (matrix [[1 1] [1 0]]) 100) (matrix* (matrix-expt (matrix [[1 1] [1 0]]) 100)
(col-matrix [0 1]))] (col-matrix [0 1]))
(->col-matrix (list (fibonacci 100) (fibonacci 99)))]
} }
@defproc[(matrix-scale [M (Matrix Number)] [z Number]) (Matrix Number)]{ @defproc[(matrix-scale [M (Matrix Number)] [z Number]) (Matrix Number)]{
@ -459,18 +476,19 @@ Returns array of the entries on the diagonal of @racket[M].
(matrix ([1 2 3] [4 5 6] [7 8 9])))] (matrix ([1 2 3] [4 5 6] [7 8 9])))]
} }
@deftogether[(@defproc[(matrix-upper-triangle [M (Matrix A)]) (Matrix (U A 0))] @deftogether[(@defproc[(matrix-upper-triangle [M (Matrix A)] [zero A 0]) (Matrix A)]
@defproc[(matrix-lower-triangle [M (Matrix A)]) (Matrix (U A 0))])]{ @defproc[(matrix-lower-triangle [M (Matrix A)] [zero A 0]) (Matrix A)])]{
The function @racket[matrix-upper-triangle] returns an upper The function @racket[matrix-upper-triangle] returns an upper
triangular matrix (entries below the diagonal are zero) with triangular matrix (entries below the diagonal have the value @racket[zero]) with
entries from the given matrix. Likewise the function entries from the given matrix. Likewise the function
@racket[matrix-lower-triangle] returns a lower triangular @racket[matrix-lower-triangle] returns a lower triangular
matrix. matrix.
@examples[#:eval untyped-eval @examples[#:eval typed-eval
(define M (array+ (array 1) (axis-index-array #(5 7) 1))) (define M (array+ (array 1) (axis-index-array #(5 7) 1)))
M M
(matrix-upper-triangle M) (matrix-upper-triangle M)
(matrix-lower-triangle M)] (matrix-lower-triangle M)
(matrix-lower-triangle (array->flarray M) 0.0)]
} }
@deftogether[(@defproc[(matrix-rows [M (Matrix A)]) (Listof (Matrix A))] @deftogether[(@defproc[(matrix-rows [M (Matrix A)]) (Listof (Matrix A))]
@ -1010,7 +1028,7 @@ The norm used by @racket[matrix-relative-error] and @racket[matrix-absolute-erro
The default value is @racket[matrix-op-inf-norm]. The default value is @racket[matrix-op-inf-norm].
Besides being a true norm, @racket[norm] should also be @deftech{submultiplicative}: Besides being a true norm, @racket[norm] should also be @deftech{submultiplicative}:
@racketblock[(norm (matrix* M0 M1)) <= (* (norm A) (norm B))] @racketblock[(norm (matrix* M0 M1)) <= (* (norm M0) (norm M1))]
This additional triangle-like inequality makes it possible to prove error bounds for formulas that This additional triangle-like inequality makes it possible to prove error bounds for formulas that
involve matrix multiplication. involve matrix multiplication.

2
collects/math/tests/matrix-strictness-tests.rkt
View File

@ -131,7 +131,7 @@
(for*: ([M (list nonstrict-2x2-arr strict-2x2-arr)]) (for*: ([M (list nonstrict-2x2-arr strict-2x2-arr)])
(check-false (array-strict? (matrix-scale M -1))) (check-false (array-strict? (matrix-scale M -1)))
(check-true (array-strict? (matrix-expt M 0))) (check-false (array-strict? (matrix-expt M 0)))
(check-false (array-strict? (matrix-expt (array-lazy M) 1))) (check-false (array-strict? (matrix-expt (array-lazy M) 1)))
(check-false (array-strict? (matrix-expt M 2))) (check-false (array-strict? (matrix-expt M 2)))
(check-false (array-strict? (matrix-expt M 3))) (check-false (array-strict? (matrix-expt M 3)))