#|

     class14.lsp   

    - Matrix factorizations for square matrices

         LU

         Cholesky

         eigen

 |#

;--- Utility for building matrix of polynomials

   (defun POLY-MATRIX (x order)
     "Builds matrix of columns x,x^2,...."
     (apply
      #'bind-columns
      (mapcar #'(lambda (p) (^ x p)) (iseq 1 order))
      ))
 
(print-matrix 
 (setf poly (poly-matrix (rseq -1 1 100) 4)))

;--- Utilities for extracting a row or column as a list

    (defun AS-LIST (x)
      (coerce (compound-data-seq x) 'list))
 
    (defun COL-AS-LIST (matrix j)
      (as-list 
       (select matrix (iseq (array-dimension matrix 0)) j)
       ))

    (defun ROW-AS-LIST (matrix i)
      (as-list 
       (select matrix i (iseq (array-dimension matrix 1)))
       ))




;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;    Decompositions: LU                                   ;;
;;                              Tierney, Section 5.5, 167  ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;--- LU factorization is most primitive, used by others (solve, detm, inv)
;    Result-> (LU in first matrix), pivoting, odd pivots, singular?

(lu-decomp '((1 2) (3 4)))       ; will not work with list

(setf lu (lu-decomp '#2a((1 2) (3 4))) )  ; nonsing

(setf lu2 (lu-decomp '#2a((0 2) (0 4))) ) ; singular


;--- Solve once, use often

(lu-solve lu '(1 6))            

(lu-solve lu '(1 3))


;--- Works with complex numbers

(setf lu (lu-decomp '#2a((#c(1 0) #c(2 0))  (#c(1 1) #c(2 3)))  ))

(lu-solve lu '(1 3))




;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;    Decompositions: Cholesky                             ;;
;;                              Tierney, Section 5.5, 168  ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;--- Most often used with psd matrices. Result is (matrix val?) list
;    Second element D is 0.0 if psd

(def cov  (covariance-matrix
           (bind-columns (normal-rand 20) (normal-rand 20))))

(setf CH (chol-decomp cov))


;--- Recover as M + D = L L'

(def c (first ch))

(print-matrix
 (- cov (matmult c (transpose c)))
 t :float-digits 2)


;--- Common application for correlated normal samples

   (defun CORR-NORMAL (n covMat)
     "Returns list of n normals with indicated covariance matrix."
     (let ((chol (chol-decomp covMat))  )
       (if (zerop (second chol))
           (let ((k (array-dimension covmat 0))  )
             (mapcar #'(lambda (x) (matmult (first chol) x))
                     (split-list (normal-rand (* n k)) k)))
           (format t "Cov not psd; offset by ~a~%" (second chol))
           )))

(def DATA (corr-normal 50 '#2a((1 .5 0) (.5 1 .5) (0 .5 1) )))

(print-matrix 
 (setf cov (apply #'covariance-matrix (transpose DATA)))
 t :float-digits 2)
     
       

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;                                                         ;;
;;    Decompositions: Eigenvalue                           ;;
;;                              Tierney, Section 5.5, 169  ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;--- Get a list of evals, evectors

(print-matrix cov)

(def e (eigen cov))

(def L (first (chol-decomp cov)))


;--- Alternative square root

(print-matrix (- cov (matmult L (transpose L))))

 
    (defun MAT-SQRT (m)
      "Eigen square root of matrix"
      (let ((e (eigen m))  )
        (apply #'+
               (mapcar #'(lambda (val vec)
                           (* (sqrt val)
                              (outer-product vec vec #'*)))
                       (coerce (first e) 'list)
                       (second e)
                       ))
        ))


(def sr (mat-sqrt cov))
(print-matrix (- cov (matmult sr (transpose sr))))