extra/math/matrices/extras/extras.factor

   1 ! Copyright (C) 2005, 2010, 2018 Slava Pestov, Joe Groff, and Cat Stevens.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: accessors arrays combinators formatting kernel math
   4 math.bits math.functions math.matrices math.order
   5 math.statistics math.text.english math.vectors random sequences
   6 sequences.deep summary ;
   7 IN: math.matrices.extras
   8
   9 ! this is a questionable implementation
  10 SINGLETONS:      +full-rank+ +half-rank+ +zero-rank+ +deficient-rank+ +uncalculated-rank+ ;
  11 UNION: rank-kind +full-rank+ +half-rank+ +zero-rank+ +deficient-rank+ +uncalculated-rank+ ;
  12
  13 ERROR: negative-power-matrix
  14     { m matrix } { n integer } ;
  15 ERROR: non-square-determinant
  16     { m integer }  { n integer } ;
  17 ERROR: undefined-inverse
  18     { m integer }  { n integer } { r rank-kind initial: +uncalculated-rank+ } ;
  19
  20 M: negative-power-matrix summary
  21     n>> dup ordinal-suffix "%s%s power of a matrix is undefined" sprintf ;
  22 M: non-square-determinant summary
  23     [ m>> ] [ n>> ] bi "non-square %s x %s matrix has no determinant" sprintf ;
  24 M: undefined-inverse summary
  25     [ m>> ] [ n>> ] [ r>> name>> ] tri "%s x %s matrix of rank %s has no inverse" sprintf ;
  26
  27 <PRIVATE
  28 DEFER: alternating-sign
  29 : finish-randomizing-matrix ( matrix -- matrix' )
  30     [ f alternating-sign randomize ] map randomize ; inline
  31 PRIVATE>
  32
  33 : <random-integer-matrix> ( m n max -- matrix )
  34     '[ _ _ 1 + random-integers ] replicate
  35     finish-randomizing-matrix ; inline
  36
  37 : <random-unit-matrix> ( m n max -- matrix )
  38     '[ _ random-units [ _ * ] map ] replicate
  39     finish-randomizing-matrix ; inline
  40
  41 <PRIVATE
  42 : (gram-schmidt) ( v seq -- newseq )
  43     [ dupd proj v- ] each ;
  44 PRIVATE>
  45
  46 : gram-schmidt ( matrix -- orthogonal )
  47     [ V{ } clone [ over (gram-schmidt) suffix! ] reduce ] keep like ;
  48
  49 : gram-schmidt-normalize ( matrix -- orthonormal )
  50     gram-schmidt [ normalize ] map ; inline
  51
  52 : kronecker-product ( m1 m2 -- m )
  53     '[ [ _ n*m  ] map ] map stitch stitch ;
  54
  55 : outer-product ( u v -- matrix )
  56     '[ _ n*v ] map ;
  57
  58 ! Special matrix constructors follow
  59 : <hankel-matrix> ( n -- matrix )
  60   [ <iota> dup ] keep '[ + abs 1 + dup _ > [ drop 0 ] when ] cartesian-map ;
  61
  62 : <hilbert-matrix> ( m n -- matrix )
  63     [ <iota> ] bi@ [ + 1 + recip ] cartesian-map ;
  64
  65 : <toeplitz-matrix> ( n -- matrix )
  66     <iota> dup [ - abs 1 + ] cartesian-map ;
  67
  68 : <box-matrix> ( r -- matrix )
  69     2 * 1 + dup '[ _ 1 <array> ] replicate ;
  70
  71 : <vandermonde-matrix> ( u n -- matrix )
  72     <iota> [ v^n ] with map reverse flip ;
  73
  74 ! Transformation matrices
  75 :: <rotation-matrix3> ( axis theta -- matrix )
  76     theta cos :> c
  77     theta sin :> s
  78     axis first3 :> ( x y z )
  79     x sq 1.0 x sq - c * +    x y * 1.0 c - * z s * -  x z * 1.0 c - * y s * + 3array
  80     x y * 1.0 c - * z s * +  y sq 1.0 y sq - c * +    y z * 1.0 c - * x s * - 3array
  81     x z * 1.0 c - * y s * -  y z * 1.0 c - * x s * +  z sq 1.0 z sq - c * +   3array
  82     3array ;
  83
  84 :: <rotation-matrix4> ( axis theta -- matrix )
  85     theta cos :> c
  86     theta sin :> s
  87     axis first3 :> ( x y z )
  88     x sq 1.0 x sq - c * +    x y * 1.0 c - * z s * -  x z * 1.0 c - * y s * +  0 4array
  89     x y * 1.0 c - * z s * +  y sq 1.0 y sq - c * +    y z * 1.0 c - * x s * -  0 4array
  90     x z * 1.0 c - * y s * -  y z * 1.0 c - * x s * +  z sq 1.0 z sq - c * +    0 4array
  91     { 0.0 0.0 0.0 1.0 } 4array ;
  92
  93 :: <translation-matrix4> ( offset -- matrix )
  94     offset first3 :> ( x y z )
  95     {
  96         { 1.0 0.0 0.0 x   }
  97         { 0.0 1.0 0.0 y   }
  98         { 0.0 0.0 1.0 z   }
  99         { 0.0 0.0 0.0 1.0 }
 100     } ;
 101
 102 <PRIVATE
 103 GENERIC: >scale-factors ( object -- x y z )
 104 M: number >scale-factors
 105     dup dup ;
 106 M: sequence >scale-factors
 107     first3 ;
 108 PRIVATE>
 109
 110 :: <scale-matrix3> ( factors -- matrix )
 111     factors >scale-factors :> ( x y z )
 112     {
 113         { x   0.0 0.0 }
 114         { 0.0 y   0.0 }
 115         { 0.0 0.0 z   }
 116     } ;
 117
 118 :: <scale-matrix4> ( factors -- matrix )
 119     factors >scale-factors :> ( x y z )
 120     {
 121         { x   0.0 0.0 0.0 }
 122         { 0.0 y   0.0 0.0 }
 123         { 0.0 0.0 z   0.0 }
 124         { 0.0 0.0 0.0 1.0 }
 125     } ;
 126
 127 : <ortho-matrix4> ( factors -- matrix )
 128     [ recip ] map <scale-matrix4> ;
 129
 130 :: <frustum-matrix4> ( xy-dim near far -- matrix )
 131     xy-dim first2 :> ( x y )
 132     near x /f :> xf
 133     near y /f :> yf
 134     near far + near far - /f :> zf
 135     2 near far * * near far - /f :> wf
 136
 137     {
 138         { xf  0.0  0.0 0.0 }
 139         { 0.0 yf   0.0 0.0 }
 140         { 0.0 0.0  zf  wf  }
 141         { 0.0 0.0 -1.0 0.0 }
 142     } ;
 143
 144 :: <skew-matrix4> ( theta -- matrix )
 145     theta tan :> zf
 146     {
 147         { 1.0 0.0 0.0 0.0 }
 148         { 0.0 1.0 0.0 0.0 }
 149         { 0.0 zf  1.0 0.0 }
 150         { 0.0 0.0 0.0 1.0 }
 151     } ;
 152
 153 ! a simpler verison of this, like matrix-map -except, but map-index, should be possible
 154 : cartesian-matrix-map ( matrix quot: ( ... pair matrix -- ... matrix' ) -- matrix-seq )
 155     [ [ first length <cartesian-square-indices> ] keep ] dip
 156     '[ _ @ ] matrix-map ; inline
 157
 158 : cartesian-column-map ( matrix quot: ( ... pair matrix -- ... matrix' ) -- matrix-seq )
 159     [ cols first2 ] prepose cartesian-matrix-map ; inline
 160
 161 ! -------------------------------------------------
 162 ! numerical analysis of matrices follows
 163 <PRIVATE
 164
 165 : square-rank ( square-matrix -- rank ) ;
 166 : nonsquare-rank ( matrix -- rank ) ;
 167 PRIVATE>
 168
 169 GENERIC: rank ( matrix -- rank )
 170 M: zero-matrix rank
 171     drop +zero-rank+ ;
 172
 173 M: square-matrix rank
 174     square-rank ;
 175
 176 M: matrix rank
 177     nonsquare-rank ;
 178
 179 GENERIC: nullity ( matrix -- nullity )
 180
 181
 182 ! implementation details of determinant and inverse
 183 <PRIVATE
 184 : alternating-sign ( seq odd-elts? -- seq' )
 185     '[ even? _ = [ neg ] unless ] map-index ;
 186
 187 ! the determinant of a 1x1 matrix is the value itself
 188 ! this works for any-dimensional matrices too
 189 : (1determinant) ( matrix -- 1det ) flatten first ; inline
 190
 191 ! optimized to find the determinant of a 2x2 matrix
 192 : (2determinant) ( matrix -- 2det )
 193     ! multiply the diagonals and subtract
 194     [ main-diagonal ] [ anti-diagonal ] bi [ first2 * ] bi@ - ; inline
 195
 196 ! optimized for 3x3
 197 ! https://www.mathsisfun.com/algebra/matrix-determinant.html
 198 :: (3determinant) ( matrix-seq -- 3det )
 199     ! first 3 elements of row 1
 200     matrix-seq first first3 :> ( a b c )
 201     ! last 2 rows, transposed to make the next step easier
 202     matrix-seq rest transpose
 203     ! get the lower sub-matrices in reverse order of a b c columns
 204     [ rest ] [ [ first ] [ third ] bi 2array ] [ 1 head* ] tri 3array
 205     ! find determinants
 206     [ (2determinant) ] map
 207     ! negate odd elements of a b c and multiply by the new determinants
 208     { a b c } t alternating-sign v*
 209     ! sum the resulting sequence
 210     sum ;
 211
 212 DEFER: (ndeterminant)
 213 : make-determinants ( n matrix -- seq )
 214     <repetition> [
 215         cols-except [ length ] keep (ndeterminant) ! recurses here
 216     ] map-index ;
 217
 218 DEFER: (determinant)
 219 ! generalized to 4 and higher
 220 : (ndeterminant) ( n matrix -- ndet )
 221     ! TODO? recurse for n < 3
 222     over 4 < [ (determinant) ] [
 223         [ nip first t alternating-sign ] [ rest make-determinants ] 2bi
 224         v* sum
 225     ] if ;
 226
 227 ! switches on dimensions only
 228 : (determinant) ( n matrix -- determinant )
 229     over {
 230         { 1 [ nip (1determinant) ] }
 231         { 2 [ nip (2determinant) ] }
 232         { 3 [ nip (3determinant) ] }
 233         [ drop (ndeterminant) ]
 234     } case ;
 235 PRIVATE>
 236
 237 GENERIC: determinant ( matrix -- determinant )
 238 M: zero-square-matrix determinant
 239     drop 0 ;
 240
 241 M: square-matrix determinant
 242     [ length ] keep (determinant) ;
 243
 244 ! determinant is undefined for m =/= n, unlike inverse
 245 M: matrix determinant
 246     dimension first2 non-square-determinant ;
 247
 248 : 1/det ( matrix -- 1/det )
 249     determinant recip ; inline
 250
 251 ! -----------------------------------------------------
 252 ! inverse operations and implementations follow
 253 ALIAS: multiplicative-inverse recip
 254
 255 ! per element, find the determinant of all other elements except the element's row / col
 256 ! https://www.mathsisfun.com/algebra/matrix-inverse-minors-cofactors-adjugate.html
 257 : >minors ( matrix -- matrix' )
 258     matrix-except-all [ [ determinant ] map ] map ;
 259
 260 ! alternately invert values of the matrix (see alternating-sign)
 261 : >cofactors ( matrix -- matrix' )
 262     [ even? alternating-sign ] map-index ;
 263
 264 ! multiply a matrix by the inverse of its determinant
 265 : m*1/det ( matrix -- matrix' )
 266     [ 1/det ] keep n*m ; inline
 267
 268 ! inverse implementation
 269 <PRIVATE
 270 ! https://www.mathsisfun.com/algebra/matrix-inverse-minors-cofactors-adjugate.html
 271 : (square-inverse) ( square-matrix -- inverted )
 272     ! inverse of the determinant of the input matrix
 273     [ 1/det ]
 274     ! adjugate of the cofactors of the matrix of minors
 275     [ >minors >cofactors transpose ]
 276     ! adjugate * 1/det
 277     bi n*m ;
 278
 279 ! TODO
 280 : (left-inverse) ( matrix -- left-invert )   ;
 281 : (right-inverse) ( matrix -- right-invert ) ;
 282
 283 ! TODO update this when rank works properly
 284 ! only defined for rank(A) = rows(A) OR rank(A) = cols(M)
 285 ! https://en.wikipedia.org/wiki/Invertible_matrix
 286 : (specialized-inverse) ( rect-matrix -- inverted )
 287     dup [ rank ] [ dimension ] bi [ = ] with map {
 288         { { t f } [ (left-inverse) ] }
 289         { { f t } [ (right-inverse) ] }
 290         [ no-case ]
 291     } case ;
 292 PRIVATE>
 293
 294 M: zero-square-matrix recip
 295     ; inline
 296
 297 M: square-matrix recip
 298     (square-inverse) ; inline
 299
 300 M: zero-matrix recip
 301     transpose ; inline ! TODO: error based on rankiness
 302
 303 M: matrix recip
 304     (specialized-inverse) ; inline
 305
 306 ! TODO: use the faster algorithm: [ determinant zero? ]
 307 : invertible-matrix? ( matrix -- ? )
 308     [ dimension first2 max <identity-matrix> ] keep
 309     dup recip mdot = ;
 310
 311 : linearly-independent-matrix? ( matrix -- ? ) ;
 312
 313 <PRIVATE
 314 ! this is the original definition of m^n as committed in 2012; it has not been lost
 315 : (m^n) ( m n -- n )
 316     make-bits over first length <identity-matrix>
 317     [ [ dupd mdot ] when [ dup mdot ] dip ] reduce nip ;
 318 PRIVATE>
 319
 320 ! A^-1 is the inverse but other negative powers are nonsense
 321 : m^n ( m n -- n ) {
 322         { [ dup -1 = ] [ drop recip ] }
 323         { [ dup 0 >= ] [ (m^n) ] }
 324         [ negative-power-matrix ]
 325     } cond ;
 326
 327 : n^m ( n m -- n ) swap m^n ; inline
 328
 329 : covariance-matrix-ddof ( matrix ddof -- cov )
 330     '[ _ cov-ddof ] cartesian-column-map ; inline
 331
 332 : covariance-matrix ( matrix -- cov )
 333     0 covariance-matrix-ddof ; inline
 334
 335 : sample-covariance-matrix ( matrix -- cov )
 336     1 covariance-matrix-ddof ; inline
 337
 338 : population-covariance-matrix ( matrix -- cov ) 0 covariance-matrix-ddof ; inline