Version 1.9 Build 1367

News	FAQ
Search	Home

Next: Ancillary Methods for Basic Data Objects Up: No Title Previous: Basic Data Values and Objects

Subsections

Linear and Non-Linear Algebra Methods and Objects

Scope

Interesting linear algebra methods begin when one considers inversions and decompositions of a matrix(A), or when one is solving linear equations like

$\displaystyle \bf A$ $\displaystyle \bf x$ = $\displaystyle \bf d$

(1)

where d is a vector of ``measurements'', and A is a ``model'' for how an vector of unknowns (x) reproduces the measurements. For non-linear models one is usually solving for a parameter set P for equations modeling measurements such as

O( $\displaystyle \sl P$ ) = $\displaystyle \bf d$

(2)

where O is some non-linear operator. Classes and methods used to determine x for linear equations, and P for non-linear equations, are what we will call solvers, and to avoid confusion we will call the generic method for solving Equation (1) linsolve and the generic method for solving Equation (2) nonlinsolve.

The algorithms for inversion, decomposition, and linsolve with a single constraint equation are reasonably well understood and we will discuss the major methods and resulting data objects in the next section. However, there are two important areas where the algorithms are not standardized for either linsolve or nonlinsolve: multiple constraint equations on the same set of unknowns; and mixtures of ``equality'' and inequality constraint equations. A simultaneous solution for unknown parameters involving four polarization equations, and inequalities imposing positivity or some other range of parameters, are cases where this would be useful. We will not say more about this here, but experimentation with multiple constraint solvers is under way. However, we can proceed knowing that all such augmentations of this type are founded on the basic data objects and the standard linear and non-linear methods that we discuss in this document.

Inversion and Decompositon

Many algorithms are based upon decomposition (also called factorization) of a square or rectangular matrices into combinations of diagonal, upper triangular, lower triangular, or other matrices. In the We propose that all decompositions produce a mlist data object with components representing the vectors or matrices produced by the decomposition, and each of these has methods specific to using that decomposition, or generic methods, used for all decompositions. Table 3 lists the most important decomposition objects with some comments about their role.

                            Table 3
              Matrix Decomposition Data Objects
Name                         Role
LUdecomp    Fastest but least robust when used in solvers,
CHOLdecomp  Choleski decomposition, best for non-zero, square
            matrices, particularly the ``normal'' equation
            matrix (tran(A) A)  
QRdecomp    Best compromise of speed and robustness when used in
            solvers, should be the default decomposition, using
            Householder transformations to perform the decomposition
SVdecomp    Singular value decomposition; slowest but most 
            robust in solvers, particularly when allowing user 
            interaction and modification of vector of singular values

Table 4 lists some of the methods that either use or support matrix decomposition objects.

                            Table 4
         Methods Related to Matrix Decomposition Data Objects
Name                         Role
backsolve     Takes an upper triangular matrix (R) and a
              vector (d) and solves for the vector (y) in R y = d
forwardsolve  Takes an lower triangular matrix (L) and a
              vector (d) and solves for the vector (x) in L x = d
invert        Takes either a decomposition object for a matrix A, 
              or A and the name of an inversion or decomposition
              method(QR, LU, Choleski, SVD), and attempts to 
              obtain Ainv such that crossprod(Ainv,tran(A)) = 1, 
              using a default or supplied tolerance value; the 
              default should be based on QR decomposition
determinant   Return a determinant of a matrix (or decomposition of
              of matrix) using a default or supplied method

The methods in Table 4 return a standard mmatrix data object and do not require creation of new data objects.

Linear Algebra Solvers and Error Analysis

The main use of matrix inversion and decomposition is in solvers with associated methods for analysis of the errors in the solution. Table 5 lists some of the data objects and methods that are basic for this form of data processing. We identify one of the matrix decomposition objects with the letter D, of any type, particularly those listed in Table 4. The letter F identifies a formula object supplied to a data fitting method.

                           Table 5
         Methods Related to Linear Solvers and Error Analysis
Name     Objects                         Role
linsolve   D        Returns the inverse of a matrix using methods
                    associated with the decomposition object
linsolve   D,V      Solves for vector x using methods of
                    the decomposition object, and the supplied 
                    ``measurement'' vector
linsolve   M        Returns the inverse of a matrix using a method
                    specified by inversion or decompostion type
linsolve   D,V      Solves for vector x using a method specified
                    by inversion or decomposition type, and the 
                    supplied ``measurement'' vector
eigen   D[,M}       returns a data object with eigenvalues and
                    eigen vectors of matrix or previously computed
                    matrix decomposition
lsfit   V,V[,V],F   Returns ``fit'' data object for given
                    (linear) formula object used in a least squares
                    fit to a ``measurement'' vector associated with
                    a vector of independent variables, and an optional
                    vector of weights, utilizing a specified 
                    decomposition method for the (tran(A) A) matrix
memfit  V,V[,V],F   Similar to lsfit except MEM is used with a
                    constraint supplied with a formula object

In Table 5 we identify a fit data object that contains the solutions to $\bf A\ x$ = $\bf D$ , and probably the residuals and the covariance matrix. The fit data objects will have additional methods return a list of traditional components of error analysis such as vectors of variance, standard deviation, correlation coefficients, etc. In the first sentence of this paragraph the word ``probably'' is used because it is unclear, when size of data sets becomes large, whether ancillary results should be automatically generated, or whether they should be generated from special methods using small fit data objects, the matrices (or decomposition objects), and the input data vectors.

There are a large number of useful fitting algorithms that could be listed in Table 5, and we are mentioning only a couple.

Please note that we are not implying that these linear algebra methods provide everything needed for AIPS++ solvers. Standard non-linear solvers, or specially coded equivalents appropriate to very large data sets are essential for many of the interesting cases; however many of these will use and operate upon basic data objects.

Number Generators

For data analysis one often needs to generate model data assuming values of the mean, standard deviation, and type of distribution (Uniform, Gaussion, Poisson, etc.). These are then used to compare real data with model distributions. Therefore part of the basic mathematics subsystem should be methods to fill VMAG data objects with model distributions of data, given specification of the desired distribution and its parameters. The basic versions already exist in AIPS++, but may need augmentation to work with the framework discussed in this document.

Statistical Analysis Methods

The statistical analysis of solutions, or fits, to data is a highly developed field for linear equations and data obeying well-understood statistics. A large number data objects and methods can be used for these standard statisical analysis methods. At first look these are the primary components of the statistics-oriented S+ . However, for the purpose of this document we eschew discussion of these approaches, assuming they will be developed as higher level applications.

Non-Linear Algebra Solvers and Error Analysis

In Table 6 we list two generic methods for solving non-linear equations. An initial version of nonlinsolve has already been protoyped in C++ using the Levenberg-Marquandt compromise for the Gauss-Newton method; however, it still needs cosmetic changes before being checked into the system. Solvers based on the Levenberg-Marquandt compromise are amongst the more robust nonlinear solvers. An early version of an iterative solver, based upon a sequence of assumed model equations supplied by formula objects, should also be prototyped early, since solvers of this type work well with most data from instruments like the VLA.

                         Table 6
              Methods Related to Non-linear Solvers

nonlinsolve   V[,M]   Solve equations of the form O(P) = d where
                      d is a supplied data vector, the form of the
                      nonlinear function O is supplied with a
                      formula object, and M is an optional matrix
                      that is part of the nonlinear function definition
  
itersolve     V[,M]   solves using iteration, and supplied formula
                      objects describing the first and subsequent
                      formulas, to solve for a parameter set

Next: Ancillary Methods for Basic Data Objects Up: No Title Previous: Basic Data Values and Objects