share/linearalgebra/linalg-extra.lisp

   1 ;; Copyright 2006 by Barton Willis
   2
   3 ;;  This is free software; you can redistribute it and/or
   4 ;;  modify it under the terms of the GNU General Public License,
   5 ;;  http://www.gnu.org/copyleft/gpl.html.
   6
   7 ;; This software has NO WARRANTY, not even the implied warranty of
   8 ;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
   9
  10 ($put '$linalgextra 1 '$version)
  11
  12 (defun $circulant (lst)
  13   ($require_list lst '$first '$circulant)
  14   (let ((q) (n ($length lst)))
  15     (setq lst (rest lst))
  16     (setq q (list lst))
  17     (decf n)
  18     (dotimes (i n)
  19       (setq lst `(,@(rest lst) ,(car lst)))
  20       (push lst q))
  21     (setq q (mapcar #'(lambda(s) (cons '(mlist) s)) q))
  22     (push '($matrix) q)))
  23
  24 (defun $cauchy_matrix (p &optional q)
  25   ($require_list p '$first '$cauchy_matrix)
  26   (if q ($require_list q '$second '$cauchy_matrix) (setq q p))
  27   (let ((row) (mat))
  28     (setq p (margs p))
  29     (setq q (margs q))
  30     (dolist (pj p)
  31       (setq row nil)
  32       (dolist (qj q)
  33         (push (div 1 (add pj qj)) row))
  34       (setq row (nreverse row))
  35       (push '(mlist) row)
  36       (push row mat))
  37     (setq mat (nreverse mat))
  38     (push '($matrix) mat)))
  39
  40 (defun $hessian (e vars)
  41   (cond (($listp vars)
  42          (let ((z) (mat nil))
  43            (setq vars (margs vars))
  44            (dolist (vi vars)
  45              (setq z ($diff e vi))
  46              (push (cons '(mlist) (mapcar #'(lambda (s) ($diff z s)) vars)) mat))
  47            (cons '($matrix) (reverse mat))))
  48         (t  `(($hessian) ,e ,vars))))
  49
  50 (defun $jacobian (e vars)
  51   (cond ((and ($listp vars) ($listp e))
  52          (setq e (margs e))
  53          (setq vars (margs vars))
  54          (let ((mat nil))
  55            (dolist (ei e)
  56              (push (cons '(mlist) (mapcar #'(lambda (s) ($diff ei s)) vars)) mat))
  57            (cons '($matrix) (reverse mat))))
  58         (t `(($jacobian) ,e ,vars))))
  59
  60 (defun $vandermonde_matrix (l)
  61   (let ((x) (q) (row) (acc nil) (n))
  62     (setq l (require-list l '$vandermonde_matrix))
  63     (setq n (- (length l) 1))
  64     (while l
  65       (setq q 1)
  66       (setq x (pop l))
  67       (setq row (list 1))
  68       (dotimes (j n)
  69         (setq q (mul q x))
  70         (push q row))
  71       (setq row (cons '(mlist) (nreverse row)))
  72       (push row acc))
  73     (simplify (cons '($matrix) (nreverse acc)))))
  74
  75 ;; Use Sylvester's criterion to decide if the self-adjoint part of a matrix is
  76 ;; negative definite (neg) or positive definite (pos). For all other cases, return
  77 ;; pnz. This algorithm is unable to determine if a matrix is negative semidefinite
  78 ;; or positive semidefinite. The argument to this function must be a selfadjoint
  79 ;; matrix.
  80
  81 (defun matrix-sign-sylvester (m)
  82   (let ((n) (det) (p-sgn nil) (n-sgn nil))
  83     (setq n ($first ($matrix_size m)))
  84
  85     (loop for i from 1 to n do
  86       (setq det (ratdisrep (newdet m i nil)))
  87       (push (csign det) p-sgn)
  88       (push (csign (mul (power -1 i) det)) n-sgn))
  89
  90     (cond ((every #'(lambda (s) (eq s '$pos)) p-sgn) '$pos)
  91           ((every #'(lambda (s) (eq s '$pos)) n-sgn) '$neg)
  92           (t '$pnz))))
  93
  94 (defun order-of-root (p x pt)
  95   (let ((order 0))
  96     (setq p ($expand p))
  97     (while (and (alike 0 ($substitute pt x p)) (not (alike 0 p)))
  98       (incf order)
  99       (setq p ($expand ($diff p x))))
 100     order))
 101
 102 ;; Let M be the self-adjoint part of mat. By the self-adjoint part
 103 ;; of a matrix M, I mean (M + M^*) / 2, where ^* is the conjugate transpose
 104 ;; Return
 105
 106 ;;  (a) neg if M is negative definite,
 107 ;;  (b) nz if M is negative semidefinite,
 108 ;;  (c) pz if M is positive semidefinite,
 109 ;;  (d) pos if M is positive definite,
 110 ;;  (e) pnz if M isn't a--d or if Maxima isn't able determine the matrix sign.
 111
 112 ;; When M is a matrix of rational numbers, look at the zeros of the characteristic polynomial;
 113 ;; otherwise, use Sylvester's criterion.
 114
 115 (defun $matrix_sign (mat)
 116   (let ((z (gensym)) (p) (n) (nbr-zero-roots) (nbr-neg-roots) (nbr-pos-roots))
 117
 118     ($require_square_matrix mat '$first '$matrix_sign)
 119     (setq mat (div (add mat ($ctranspose mat)) 2))
 120     (cond (($some '((lambda) ((mlist) $s) ((mnot) (($ratnump) $s))) mat)
 121            (matrix-sign-sylvester mat))
 122           (t
 123            (setq n ($first ($matrix_size mat)))
 124            (cond (($zeromatrixp mat) '$zero)
 125                  (t
 126                   (setq p ($charpoly mat z))
 127
 128                   ;; number of roots of characteristic poly in {0}
 129                   (setq nbr-zero-roots (order-of-root p z 0))
 130
 131                   ;; number of roots of characteristic poly in (-inf,0)
 132                   (setq nbr-neg-roots (- ($nroots p '$minf 0) nbr-zero-roots))
 133
 134                   ;; number of roots of characteristic poly in (0,inf)
 135                   (setq nbr-pos-roots ($nroots p 0 '$inf))
 136
 137                   (cond ((= n nbr-neg-roots) '$neg)
 138                         ((= 0 nbr-pos-roots) '$nz)
 139                         ((= n nbr-pos-roots) '$pos)
 140                         ((= n (+ nbr-pos-roots nbr-zero-roots)) '$pz)
 141                         (t '$pnz))))))))
 142
 143 (defun $sylvester_matrix (p q z)
 144   (let ((p-coeff nil) (q-coeff nil) (mat nil) (p-deg) (q-deg))
 145     (if (or (not ($polynomialp p `((mlist) ,z) `((lambda) ((mlist) c) (($freeof) ,z c))))
 146             (not ($polynomialp q `((mlist) ,z) `((lambda) ((mlist) c) (($freeof) ,z c)))))
 147         (merror "The first two arguments to 'sylvester_matrix' must be polynomials"))
 148
 149     (setq p ($ratexpand p))
 150     (setq q ($ratexpand q))
 151
 152     (setq p-deg ($hipow p z))
 153     (setq q-deg ($hipow q z))
 154
 155     (dotimes (k (- q-deg 1))
 156       (push 0 p-coeff))
 157
 158     (dotimes (k (+ p-deg 1))
 159       (push ($ratcoef p z k) p-coeff))
 160
 161     (dotimes (k (- p-deg 1))
 162       (push 0 q-coeff))
 163
 164     (dotimes (k (+ 1 q-deg))
 165       (push ($ratcoef q z k) q-coeff))
 166
 167     (dotimes (k q-deg)
 168       (push (cons '(mlist) p-coeff) mat)
 169       (setq p-coeff (butlast p-coeff))
 170       (push 0 p-coeff))
 171
 172     (dotimes (k p-deg)
 173       (push (cons '(mlist) q-coeff) mat)
 174       (setq q-coeff (butlast q-coeff))
 175       (push 0 q-coeff))
 176
 177     (simplify `(($matrix) ,@(reverse mat)))))
 178
 179 (defun $krylov_matrix (v mat &optional (n 'no-value))
 180   ($require_square_matrix mat '$second '$krylov_matrix)
 181   (and ($listp v) (setq v ($transpose v)))
 182   (if (not
 183        (and ($matrixp v)
 184        (= ($first ($matrix_size v)) ($second ($matrix_size mat)))))
 185       (merror "Incompatible matrix sizes"))
 186   (if (eq n 'no-value) (setq n ($first ($matrix_size mat))))
 187   ($require_posinteger n  '$third '$krylov_matrix)
 188   (let ((acc v))
 189     (decf n 1)
 190     (dotimes (k n)
 191       (setq v (ncmul mat v))
 192       (setq acc ($addcol acc v)))
 193     acc))