lsfit - Tool

Postscript file available

Least-squares fit by matrix inversion

include "lsfit.g"

Constructors
Functions

accumulate	accumulate an equation
help	print summary of lsfit functions
init	initialise before accumulation
inityonly	initialise the vector only
print	print internal matrices
solve	solves for the unknowns (x)
status	print lsfit internal status
test	test of the main lsfit functions

Description

The functions in this tool have not been tested as rigorously as other parts of AIPS++. This tool can be considered ``alpha'' code.

lsfit The main function of this tool is simultaneous solution of a set of linear equations by matrix inversion. If the matrix is rectangular, i.e. if there are more equations than unknowns, this amounts to a least-squares fit.

$\displaystyle \vec{y}\,$	=	A $\displaystyle \vec{x}\,$	(1.22)
A^*T $\displaystyle \vec{y}\,$	=	(A^*TA) $\displaystyle \vec{x}\,$	(1.23)

in which (*) indicates complex conjugation. The latter expression represents a set of normal equations, which can be solved by matrix inversion:

$\displaystyle \vec{x}\,$ = (A^*TA)^-1 (A^*T $\displaystyle \vec{y}\,$ )

(1.24)

Various matrix inversion functions are supported: Gauss-Jordan, Choleski and Singular Value Decomposition (SVD). The latter automatically takes care of singular matrices, i.e. those cases where one is trying to solve for unknowns without supplying sufficient information. The default inversion function is Choleski decomposition, because of its efficiency. The procedure is as follows (see also the example below): It is assumed that the number of unknowns (i.e. the length of the vector $\vec{x}\,$ ) is known, but that the equations $\vec{y}\,$ = A $\vec{x}\,$ will become available one by one. First, an n x n matrix (A^*TA) and a vector (A^*T $\vec{y}\,$ ) with length n are initialised with zeroes. Then, as each equation i becomes available, its coefficients a_ij and `driving' value y_i are accumulated in the matrix and the vector. The system can be solved by matrix inversion as soon as at least n equations have been accumulated.

Since the solution does not destroy the accumulation matrix and vector, one may continue to accumulate more equations after a solve. In this way, one may get intermediate results.

If one needs more than one solution for different driving vectors $\vec{y}\,$ , but the same matrix coefficients a_ij, one may use the same inverted matrix each time, and only re-accumulate the vector (A^*T $\vec{y}\,$ ). This saves time if the n is large.

The lsfit tool is a place-holder, until this functionality will be provided by a fitter DO. In the meantime, it is a good example of the power of prototyping in Glish, to find out what functionality and interface is really needed.

Example

NB: In view of the loop, it is recommended to make a .g script to try this example.

include "lsfit.g"
nuk := 5;                         # nr of unknowns		
print 'xxin=',xxin := [1:nuk];    # simulated input values 
lsf.init(nuk);	                  # initialise matrix and vector
print lsf.status();               # optional
for (i in [1:(nuk-1)]) {
    for (j in [(i+1):nuk]) {
        aa := rep(0,nuk);         # coefficient vector
        aa[i] := 1;
        aa[j] := 1;	
        y := xxin[i]*aa[i] + xxin[j]*aa[j]; # driving value
        lsf.accumulate(aa,y);     # accumulate equation
    }
}
print 'xxout=',xxout := lsf.solve(); # estimated values

Next: lsfit.init - Function Up: mathematics - Module Previous: matrix.test - Function Contents Index

	Package	utility
	Module	mathematics