newmata.txt

矩阵计算库

        ====  == =============

This is a variant on the usual QR transformation.

Start with matrix

       / 0    0 \      s
       \ X    Y /      n

         s    t

Our version of the QR decomposition multiplies this matrix by an orthogonal
matrix Q to get

       / U    M \      s
       \ 0    Z /      n

         s    t

where 'U' is upper triangular (the R of the QR transform).

This is good for solving least squares problems: choose b (matrix or row
vector) to minimize the sum of the squares of the elements of

         Y - X*b

Then choose 'b = U.i()*M;' The residuals 'Y - X*b' are in 'Z'.

This is the usual QR transformation applied to the matrix 'X' with the square
zero matrix attached concatenated on top of it. It gives the same triangular
matrix as the QR transform applied directly to 'X' and generally seems to work
in the same way as the usual QR transform. However it fits into the matrix
package better and also gives us the residuals directly. It turns out to be
essentially a modified Gram-Schmidt decomposition.

Two routines are provided:

    QRZ(X, U);

replaces 'X' by orthogonal columns and forms 'U'.

    QRZ(X, Y, M);

uses 'X' from the first routine, replaces 'Y' by 'Z' and forms 'M'.

The are also two routines 'QRZT(X, L)' and 'QRZT(X, Y, M)' which do the same
decomposition on the transposes of all these matrices. QRZT replaces the
routines HHDecompose in earlier versions of newmat. HHDecompose is still
defined but just calls QRZT.



        3.19  Singular value decomposition
        ====  ======== ===== =============

The singular value decomposition of an m x n matrix 'A' (where m >= n) is a
decomposition

    A  = U * D * V.t()

where 'U' is m x n with  'U.t() * U'  equalling the identity, 'D' is an n x n
diagonal matrix and 'V' is an n x n orthogonal matrix.

Singular value decompositions are useful for understanding the structure of
ill-conditioned matrices, solving least squares problems, and for finding the
eigenvalues of 'A.t() * A'.

To calculate the singular value decomposition of 'A' (with m >= n) use one of

    SVD(A, D, U, V);                  // U (= A is OK)
    SVD(A, D);
    SVD(A, D, U);                     // U (= A is OK)
    SVD(A, D, U, FALSE);              // U (can = A) for workspace only
    SVD(A, D, U, V, FALSE);           // U (can = A) for workspace only

The values of 'A' are not changed unless 'A' is also inserted as the third
argument.



        3.20  Eigenvalue decomposition
        ====  ========== =============

An eigenvalue decomposition of a symmetric matrix 'A' is a decomposition

    A  = V * D * V.t()

where 'V' is an orthogonal matrix and 'D' is a diagonal matrix.

Eigenvalue analyses are used in a wide variety of engineering, statistical and
other mathematical analyses.

The package includes two algorithms: Jacobi and Householder. The first is
extremely reliable but much slower than the second.

The code is adapted from routines in "Handbook for Automatic Computation, Vol
II, Linear Algebra" by Wilkinson and Reinsch, published by Springer Verlag. 


    Jacobi(A,D,S,V);                  // A, S symmetric; S is workspace,
                                      //    S = A is OK; V is a matrix
    Jacobi(A,D);                      // A symmetric
    Jacobi(A,D,S);                    // A, S symmetric; S is workspace,
                                      //    S = A is OK
    Jacobi(A,D,V);                    // A symmetric; V is a matrix

    EigenValues(A,D);                 // A symmetric
    EigenValues(A,D,S);               // A, S symmetric; S is for back
                                      //    transforming, S = A is OK
    EigenValues(A,D,V);               // A symmetric; V is a matrix

The values of 'A' are not changed unless 'A' is also inserted as the third
argument. If you need eigenvectors use one of the forms with matrix 'V'. The
eigenvectors are returned as the columns of 'V'.




        3.21  Sorting
        ====  =======

To sort the values in a matrix or vector, 'A', (in general this operation
makes sense only for vectors and diagonal matrices) use

    SortAscending(A);

or

    SortDescending(A);

I use the Shell-sort algorithm. This is a medium speed algorithm, you might
want to replace it with something faster if speed is critical and your
matrices are large.



        3.22  Fast Fourier transform
        ====  ==== ======= =========


   FFT(X, Y, F, G);                         // F=X and G=Y are OK

where X, Y, F, G are column vectors. X and Y are the real and imaginary input
vectors; F and G are the real and imaginary output vectors. The lengths of X
and Y must be equal and should be the product of numbers less than about 10
for fast execution.

The formula is

          n-1
   h[k] = SUM  z[j] exp (-2 pi i jk/n)
          j=0

where z[j] is stored complex and stored in X(j+1) and Y(j+1). Likewise h[k] is
complex and stored in F(k+1) and G(k+1). The fast Fourier algorithm takes
order n log(n) operations (for "good" values of n) rather than n**2 that
straight evaluation would take.

I use the method of Carl de Boor (1980), "Siam J Sci Stat Comput", pp 173-8.
The sines and cosines are calculated explicitly. This gives better accuracy,
at an expense of being a little slower than is otherwise possible.

Related functions

   FFTI(F, G, X, Y);                        // X=F and Y=G are OK
   RealFFT(X, F, G);
   RealFFTI(F, G, X);

'FFTI' is the inverse trasform for 'FFT'. 'RealFFT' is for the case when the
input vector is real, that is Y = 0. I assume the length of X, denoted by n,
is even. The program sets the lengths of F and G to n/2 + 1. 'RealFFTI' is the
inverse of 'RealFFT'.



        3.23  Interface to Numerical Recipes in C
        ====  ========= == ========= ======= == =

This package can be used with the vectors and matrices defined in "Numerical
Recipes in C". You need to edit the routines in Numerical Recipes so that the
elements are of the same type as used in this package. Eg replace float by
double, vector by dvector and matrix by dmatrix, etc. You will also need to
edit the function definitions to use the version acceptable to your compiler.
Then enclose the code from Numerical Recipes in  extern "C" { ... }. You will
also need to include the matrix and vector utility routines.

Then any vector in Numerical Recipes with subscripts starting from 1 in a
function call can be accessed by a RowVector, ColumnVector or DiagonalMatrix
in the present package. Similarly any matrix with subscripts starting from 1
can be accessed by an  nricMatrix  in the present package. The class
nricMatrix is derived from Matrix and can be used in place of Matrix. In each
case, if you wish to refer to a RowVector, ColumnVector, DiagonalMatrix or
nricMatrix X in an function from Numerical Recipes, use X.nric() in the
function call.

Numerical Recipes cannot change the dimensions of a matrix or vector. So
matrices or vectors must be correctly dimensioned before a Numerical Recipes
routine is called.

For example

   SymmetricMatrix B(44);
   .....                             // load values into B
   nricMatrix BX = B;                // copy values to an nricMatrix
   DiagonalMatrix D(44);             // Matrices for output
   nricMatrix V(44,44);              //    correctly dimensioned
   int nrot;
   jacobi(BX.nric(),44,D.nric(),V.nric(),&nrot);
                                     // jacobi from NRIC
   cout << D;                        // print eigenvalues




        3.24  Exceptions
        ====  ==========

This package includes a partial implementation of exceptions for compilers
which do not support C++ exceptions. I used Carlos Vidal's article in the
September 1992 "C Users Journal" as a starting point.

See customising [2.2] for selecting the correct exception option.

See the section on error messages [4] for some notes on the messages generated
by the exceptions.

The rest of this section describes my exception simulation package. But see
the end of the section for controlling the action after an exception has been
generated.

Newmat does a partial clean up of memory following throwing an exception - see
the next section. However, the present version will leave a little heap memory
unrecovered under some circumstances. I would not expect this to be a major
problem, but it is something that needs to be sorted out.

The functions/macros I define are Try, Throw, Catch, CatchAll and
CatchAndThrow. Try, Throw, Catch and CatchAll correspond to try, throw, catch
and catch(...) in the C++ standard. A list of Catch clauses must be terminated
by either CatchAll or CatchAndThrow but not both. Throw takes an Exception as
an argument or takes no argument (for passing on an exception). I do not have
a version of Throw for specifying which exceptions a function might throw.
Catch takes an exception class name as an argument; CatchAll and CatchAndThrow
don't have any arguments. Try, Catch and CatchAll must be followed by blocks
enclosed in curly brackets.

I have added another macro Bounce to mean a rethrow, Throw(). This was
necessary to enable the package to be compatible with both my exception
package and C++ exceptions.

All Exceptions must be derived from a class, Exception, defined in newmat and
can contain only static variables. See the examples in newmat if you want to
define additional exceptions.

I have defined 5 clases of exceptions for users (there are others but I
suggest you stick to these ones):

   SpaceException                 Insufficient space on the heap
   ProgramException               Errors such as out of range index or
                                  incompatible matrix types or
                                  dimensions
   ConvergenceException           Iterative process does not converge
   DataException                  Errors such as attempting to invert a
                                  singular matrix
   InternalException              Probably a programming error in newmat

For each of these exception classes, I have defined a member function 'void
SetAction(int)'. If you call 'SetAction(1)', and a corresponding exception
occurs, you will get an error message. If there is a Catch clause for that
exception, execution will be passed to that clause, otherwise the program will
exit. If you call 'SetAction(0)' you will get the same response, except that
there will be no error message. If you call 'SetAction(-1)', you will get the
error message but the program will always exit.



        3.25  Cleanup after an exception
        ====  ======= ===== == =========

The exception mechanisms in newmat are based on the C functions setjmp and
longjmp. These functions do not call destructors so can lead to garbage being
left on the heap. (I refer to memory allocated by "new" as heap memory). For
example, when you call

   Matrix A(20,30);

a small amount of space is used on the stack containing the row and column
dimensions of the matrix and 600 doubles are allocated on the heap for the
actual values of the matrix. At the end of the block in which A is declared,
the destructor for A is called and the 600 doubles are freed. The locations on
the stack are freed as part of the normal operations of the stack. If you
leave the block using a longjmp command those 600 doubles will not be freed
and will occupy space until the program terminates.

To overcome this problem newmat keeps a list of all the currently declared
matrices and its exception mechanism will return heap memory when you do a
Throw and Catch.

However it will not return heap memory from objects from other packages. If
you want the mechanism to work with another class you will have to do four
things:

1:  derive your class from class Janitor defined in except.h;

2:  define a function 'void CleanUp()' in that class to return all heap
memory;

3:   include the following lines in the class definition

      public:
         void* operator new(size_t size)
         { do_not_link=TRUE; void* t = ::operator new(size); return t; }
         void operator delete(void* t) { ::operator delete(t); }

4:   be sure to include a copy constructor in you class definition, that is,
something like

      X(const X&);


Note that the function 'CleanUp()' does somewhat the same duties as the
destructor. However 'CleanUp()' has to do the "cleaning" for the class you are
working with and also the classes it is derived from. So it will often be
wrong to use exactly the same code for both 'CleanUp()' and the destructor or
to define your destructor as a call to 'CleanUp()'.




        3.26  Non-linear applications
        ====  ========== ============

Files solution.h, solution.cpp contain a class for solving for x in y = f(x)
where x is a one-dimensional continuous monotonic function. This is not a
matrix thing at all but is included because it is a useful thing and because
it is a simpler version of the technique used in the non-linear least squares.

Files newmatnl.h, newmatnl.cpp contain a series of classes for non-linear
least squares, to be extended to maximum likelihood in a later release.

Documentation for both of these is in the definition files. Simple examples
are in sl_ex.cpp and nl_ex.cpp.

These packages do not work with the AT&T [2.3.1] and HPUX [2.3.4] compilers
and newmatnl is questionabl