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