(FILECREATED "29-May-86 12:41:55" {QV}<PEDERSEN>LISP>OTHERMATRIXOPS.;2 19199
changes to: (VARS OTHERMATRIXOPSCOMS)
(FNS SVDNASH)
previous date: "27-May-86 22:00:09" {QV}<PEDERSEN>LISP>OTHERMATRIXOPS.;1)
(* Copyright (c) 1986 by Xerox Corporation. All rights reserved.)
(PRETTYCOMPRINT OTHERMATRIXOPSCOMS)
(RPAQQ OTHERMATRIXOPSCOMS ((* File created by Coms Manager.)
(FNS DIAGXTXINVERSE OLSMGS QRMGS SVDNASH)))
(* File created by Coms Manager.)
(DEFINEQ
(DIAGXTXINVERSE
[LAMBDA (R D) (* jop: " 3-Feb-86 15:41")
(* * Return (as the P vector D) the diagonal elements of the matrix
(R-transpose R) Inverse, where R is P x P upper triangular with non-zero
diagonal elements. When the decomposition X=Q*R has been formed, these diagonal
elements are just those of (X-transpose X) Inverse
(mathematically)%. Notice that the vector D must be created by the CALLER!)
(BLAS.CHECKARRAY R)
(BLAS.CHECKARRAY D)
(PROG [(P (ARRAY-DIMENSION R 0))
(C (MAKE-ARRAY (ARRAY-DIMENSION R 0)
(QUOTE :ELEMENT-TYPE)
(QUOTE FLOAT]
(* * Check args)
(if (NOT (IEQP P (ARRAY-DIMENSION R 1)))
then (HELP "R not P x P" R))
(if (NOT (IEQP P (ARRAY-DIMENSION D 0)))
then (HELP "D not a P vector" D))
(* * The solution proceeds as a succession of back substitutions against
columns of the identity matrix.)
(BLAS.ARRAYFILL 0.0 D)
(bind FTEMP declare (TYPE FLOATP FTEMP) for K from (SUB1 P) to 0 by -1
do (\FLOATASET (FQUOTIENT 1.0 (\FLOATAREF R K K))
C K)
(SETQ FTEMP (\FLOATAREF C K))
(\FLOATASET (FPLUS (\FLOATAREF D K)
(FTIMES FTEMP FTEMP))
D K)
(for J from (SUB1 K) to 0 by -1
do (SETQ FTEMP 0.0)
(for I from (ADD1 J) to K do [SETQ FTEMP (FDIFFERENCE FTEMP
(FTIMES (\FLOATAREF R J I)
(\FLOATAREF C I]
finally (\FLOATASET (FQUOTIENT FTEMP (\FLOATAREF R J J))
C J))
(SETQ FTEMP (\FLOATAREF C J))
(\FLOATASET (FPLUS (\FLOATAREF D J)
(FTIMES FTEMP FTEMP))
D J)))
(RETURN NIL])
(OLSMGS
[LAMBDA (Y X Q R B) (* jop: "27-May-86 17:08")
(* * OLSMGS calculates the least squares
(multiple) regression of Y on X. An N vector Y, and a N x P matrix X are
required. Q is a N x P matrix containing an orthogonal basis for the column
spce of X. R is an upper triangular P x P array.
B is a P vector of coefficients)
(PROG ((N (ARRAY-DIMENSION X 0))
(P (ARRAY-DIMENSION X 1))
C)
(* * Arg checks)
(if (NOT (IEQP (ARRAY-TOTAL-SIZE Y)
N))
then (HELP "Y not an N vector" Y))
(if (NULL Q)
then (SETQ Q (MAKE-ARRAY (ARRAY-DIMENSIONS X)
(QUOTE :ELEMENT-TYPE)
(QUOTE FLOAT)))
elseif (NOT (EQUAL (ARRAY-DIMENSIONS Q)
(ARRAY-DIMENSIONS X)))
then (HELP "Q not N x P" Q))
(if (NULL R)
then (SETQ R (MAKE-ARRAY (LIST P P)
(QUOTE :ELEMENT-TYPE)
(QUOTE FLOAT)))
elseif (NOT (EQUAL (ARRAY-DIMENSIONS R)
(LIST P P)))
then (HELP "R not P x P" R))
(if (NULL B)
then (SETQ B (MAKE-ARRAY P (QUOTE :ELEMENT-TYPE)
(QUOTE FLOAT)))
elseif (NOT (IEQP (ARRAY-TOTAL-SIZE R)
P))
then (HELP "B not a P vector" B))
(SETQ C (MAKE-ARRAY P (QUOTE :ELEMENT-TYPE)
(QUOTE FLOAT)))
(* * Compute Q R decomposition)
(QRMGS X Q R)
(for J from 0 to (SUB1 P) do (bind (FTEMP ← 0.0) declare (TYPE FLOATP FTEMP) for I
from 0 to (SUB1 N)
do [SETQ FTEMP (FPLUS FTEMP (FTIMES (\FLOATAREF Q I J)
(\FLOATAREF Y I]
finally (\FLOATASET FTEMP C J)))
(RSOLV R C B)
(RETURN B])
(QRMGS
[LAMBDA (X Q R) (* jop: " 3-Feb-86 16:04")
(* * This function determines a decomposition X = Q*R for the rectangular
CMLArray matrix X, where Q is an orthonormal basis for the column space of X
and R is a square upper triangular matrix.)
(* * On entry: X is an n x p CMLARRAY, i.e.
a matrix, whose QR decomposition is required;
Q is an n x p array; R is a p x p array)
(* * On return: the value of the function is NIL for a normal return.
Should this routine be unable to complete the QR decomposition of X because of
near singularity of its columns, the function will have a nonNIL value;
the matrix Q has been OVERWRITTEN by the orthonormal basis Q;
R contains in its upper triangle the factor R from the decomposition.)
(* * Algorithm notes: The method employed for the decomposition is
modified-Gram-Schmidt with reorthogonalization of the columns of Q.
A test relative to the machine arithmetic precision is made to recognize
possible singularity of X. The decision to reorthoganalize is based on the
relative decrease in column norm obtained in the Gram-Schmidt process.
These ideas are due to David M. Gay. For the sake of time and space economy,
this function does not perform column pivoting;
this will exclude successful decomposition for only the most severely
ill-conditioned problems. For the sake of time and space economy, no special
precautions have been taken in the computation of inner products and vector
2-norms. Only those problems where intermediate quantities lie outside the
range 10↑-19 to 10↑19 will notice difficulties, but values of these magnitudes
are surely indication of very poor scaling or severe ill-conditioning.)
(BLAS.CHECKARRAY X)
(BLAS.CHECKARRAY Q)
(BLAS.CHECKARRAY R)
(PROG [(N (ARRAY-DIMENSION X 0))
(P (ARRAY-DIMENSION X 1))
(C (MAKE-ARRAY (ARRAY-DIMENSION X 1)
(QUOTE :ELEMENT-TYPE)
(QUOTE FLOAT]
(* * Check args)
(if (NOT (AND (IEQP (ARRAY-DIMENSION Q 0)
N)
(IEQP (ARRAY-DIMENSION Q 1)
P)))
then (HELP "Q not n x p" Q))
(if (NOT (AND (IEQP (ARRAY-DIMENSION R 0)
P)
(IEQP (ARRAY-DIMENSION R 1)
P)))
then (HELP "R not p x p" R))
(* * Copy X to Q)
(BLAS.ARRAYBLT X 0 1 Q 0 1)
(* * Initialize diagonal of R to 2-norms of columns of Q)
(bind NORMXJ for J from 0 to (SUB1 P) do (SETQ NORMXJ (BLAS.NRM2 Q J P N))
(\FLOATASET NORMXJ R J J)
(\FLOATASET (if (FGREATERP NORMXJ 0.0)
then 0.0
else 1.0)
C J))
(* Loop over each column of Q in turn.)
[bind RKK CK V8 V9 declare (TYPE FLOATP RKK CK V8) for K from 0 to (SUB1 P)
do (SETQ RKK (\FLOATAREF R K K))
(\FLOATASET 1.0 R K K)
(SETQ V9 RKK)
(SETQ CK (\FLOATAREF C K))
(bind NO.ORTHOFLG repeatuntil NO.ORTHOFLG
do [if (UFGREATERP CK .75)
then (* Weight loss has been substantial,
recompute 2-norm.)
(SETQ V8 (BLAS.NRM2 Q K P N))
else (* Otherwise, fast update of length.)
(SETQ V8 (FTIMES V9 (SQRT (SETQ V8 (FDIFFERENCE 1.0 CK]
(\FLOATASET (FTIMES V8 (\FLOATAREF R K K))
R K K)
(* Test if new column length is so much smaller than original length that we
must give up and declare singularity. The test relies implicitly on the finite
precision of machine addition.)
(if [NOT (UFLESSP RKK (FPLUS RKK (\FLOATAREF R K K]
then (RETURN K)
else (* Otherwise, accept column length.)
(for I from 0 to (SUB1 N)
do (\FLOATASET (FQUOTIENT (\FLOATAREF Q I K)
V8)
Q I K))) (* Column K should be nearly
orthogonal to the previous ones, and
we check if it can benefit from
reorthoganalization.)
(if (UFLESSP (FQUOTIENT V8 V9)
.25)
then (SETQ NO.ORTHOFLG NIL) (* Reorthoganalize with respect
preceding columns.)
(bind DOTXJXK declare (TYPE FLOATP DOTXJXK) for J from 0
to (SUB1 K)
do (SETQ DOTXJXK (BLAS.DOTPROD Q J P Q K P N))
(\FLOATASET (FPLUS (\FLOATAREF R J K)
(FTIMES DOTXJXK (\FLOATAREF R K K)))
R J K)
(for I from 0 to (SUB1 N)
do (\FLOATASET (FDIFFERENCE (\FLOATAREF Q I K)
(FTIMES DOTXJXK (\FLOATAREF Q I J)))
Q I K)))
(SETQ V9 1.0)
else (SETQ NO.ORTHOFLG T))) (* Orthogonalize remaining columns
with respect to current one.)
(bind DOTXJXK declare (TYPE FLOATP DOTXJXK) for J from (ADD1 K)
to (SUB1 P) do (SETQ DOTXJXK (BLAS.DOTPROD Q J P Q K P N))
(\FLOATASET DOTXJXK R K J)
(for I from 0 to (SUB1 N)
do (\FLOATASET (FDIFFERENCE (\FLOATAREF Q I J)
(FTIMES DOTXJXK (\FLOATAREF Q I K)))
Q I J))
(SETQ V9 (\FLOATAREF R J J))
(if (UFGREATERP V9 0.0)
then (\FLOATASET (FPLUS (\FLOATAREF C J)
(FTIMES (FQUOTIENT DOTXJXK V9)
(FQUOTIENT DOTXJXK V9)))
C J]
(RETURN NIL])
(SVDNASH
[LAMBDA (UMATRIX SVECTOR VMATRIX) (* jop: "31-Jan-86 13:50")
(* * Singular-value decomposition by means of orthogonalization by plane
rotations. Taken from Nash and Shlien: "Partial svd algorithms." On entry
UMATRIX contains the m by n matrix to be decomposed, SVECTOR must be a vector
of length n, and VMATRIX must be a square n by n matrix.
On return UMATRIX has been overwritten by the left singular vectors, SVECTOR
contains the singular values, and VMATRIX contains the right singular vectors.)
(LET* ((M (ARRAY-DIMENSION UMATRIX 0))
(N (ARRAY-DIMENSION UMATRIX 1))
(EPS 1.0E-6)
(SLIMIT (IMAX (IQUOTIENT N 4)
6))
(NT N))
(* * Initialize VMATRIX to identity matrix.)
(BLAS.ARRAYFILL 0.0 VMATRIX)
(for I from 0 to (SUB1 N) do (LASET 1.0 VMATRIX I I))
(* * The main loop: repeatedly sweep over all pairs of columns in U, rotating
as needed, until no rotations in a complete sweep are effective.
Check the opportunity for rank reduction at the conclusion of each sweep.)
[bind (SCOUNT ← 0)
RCOUNT
(FEPS ← EPS)
(FEPSSQUARED ← (FTIMES EPS EPS)) declare (TYPE FLOATP FEPS FEPSSQUARED)
repeatwhile (IGREATERP RCOUNT 0) eachtime (SETQ RCOUNT (IQUOTIENT (ITIMES NT
(SUB1 NT))
2))
(SETQ SCOUNT (ADD1 SCOUNT))
do (if (IGREATERP SCOUNT SLIMIT)
then (HELP "Number of sweeps exceeds sweep limit." SCOUNT))
[for J from 0 to (IDIFFERENCE NT 2)
do (bind P Q R C FC S FS V for K from (IPLUS J 1) to (SUB1 NT)
declare (TYPE FLOATP P Q R FC FS V)
do (SETQ P (BLAS.DOTPROD UMATRIX J N UMATRIX K N M))
(SETQ Q (BLAS.DOTPROD UMATRIX J N UMATRIX J N M))
(SETQ R (BLAS.DOTPROD UMATRIX K N UMATRIX K N M))
(LASET Q SVECTOR J)
(LASET R SVECTOR K)
(if (UFLESSP Q R)
then (SETQ Q (FDIFFERENCE (FQUOTIENT Q R)
1.0))
(SETQ P (FQUOTIENT P R))
[SETQ V (SQRT (SETQ V (FPLUS (FTIMES 4.0 P P)
(FTIMES Q Q]
[SETQ FS (SQRT (SETQ FS (FTIMES .5 (FDIFFERENCE 1.0
(FQUOTIENT Q V]
(if (UFLESSP P 0.0)
then (SETQ FS (FDIFFERENCE 0.0 FS)))
(SETQ FC (FQUOTIENT P (FTIMES V FS)))
(* box before the COLROT calls)
(SETQ C FC)
(SETQ S FS)
(BLAS.ROT C S UMATRIX J N UMATRIX K N M)
(BLAS.ROT C S VMATRIX J N VMATRIX K N N)
elseif (OR (UFLEQ (FTIMES Q R)
FEPSSQUARED)
(UFLEQ (FTIMES (FQUOTIENT P Q)
(FQUOTIENT P R))
FEPS))
then (SETQ RCOUNT (SUB1 RCOUNT))
else (SETQ R (FDIFFERENCE 1.0 (FQUOTIENT R Q)))
(SETQ P (FQUOTIENT P Q))
[SETQ V (SQRT (SETQ V (FPLUS (FTIMES 4.0 P P)
(FTIMES R R]
[SETQ FC (SQRT (SETQ FC (FTIMES .5 (FPLUS 1.0 (FQUOTIENT R V]
(SETQ FS (FQUOTIENT P (FTIMES V FC)))
(* box before the COLROT calls)
(SETQ C FC)
(SETQ S FS)
(BLAS.ROT C S UMATRIX J N UMATRIX K N M)
(BLAS.ROT C S VMATRIX J N VMATRIX K N N]
(bind (TOL ← EPS) declare (TYPE FLOATP TOL)
while (AND (IGEQ NT 3)
(UFLEQ (FQUOTIENT (LAREF SVECTOR (SUB1 NT))
(FPLUS (LAREF SVECTOR 0)
TOL))
TOL)) do (SETQ NT (SUB1 NT]
(* * Finish the decomposition by returning all N singular values, and by
normalizing those columns of UMATRIX judged to be non-zero.)
(bind Q for J from 0 to (SUB1 N) do (SETQ Q (SQRT (LAREF SVECTOR J)))
(LASET Q SVECTOR J)
(if (ILEQ J NT)
then (BLAS.SCAL (FQUOTIENT 1.0 Q)
UMATRIX J N M])
)
(PUTPROPS OTHERMATRIXOPS COPYRIGHT ("Xerox Corporation" 1986))
(DECLARE: DONTCOPY
(FILEMAP (NIL (529 19114 (DIAGXTXINVERSE 539 . 2873) (OLSMGS 2875 . 5265) (QRMGS 5267 . 13239) (
SVDNASH 13241 . 19112)))))
STOP