SVDDIV

Purpose:

Performs general least squares fitting of basis functions using singular value decomposition.

Syntax:

SVDDIV(A, y, tol)

(coef, coefstd) = SVDDIV(A, y, tol)

A

-

An input array, the basis function design matrix.

y

-

An input series, the Y output.

tol

-

Optional. A real, the threshold at which to zero out singular values. If not specified, all singular values are used.

Returns:

A series, the fitted coefficients.

 

(coef, coefstd) = SVDDIV(A, y, tol) returns the fitted coefficients and the standard deviation of the coefficients.

Example:

x = {-1, 1, 1.5, 3};

y = {6, 2, 2.25, 6};

 

c1 = polyfit(y, x, 2);

c2 = svddiv(ravel(ones(length(x),1), x, x^2), y);

 

c1 == c2 == {3, -2, 1};

 

The first method uses POLYFIT to fit a quadratic polynomial to the data points. As shown, SVDDIV produces the same result. The fitted equation becomes:

 

image\svddiv27.gif

Example:

W1: gnorm(100,.01)

W2: 5*W1 + 3*sin(W1) - 2*exp(W1*W1)

W3: ravel(W1, sin(W1), exp(W1*W1))

W4: svddiv(W3, W2)

 

W4 == {5, 3, -2}, the coefficients of W2.

 

In this example, we wish to find the coefficients of the relation:

 

image\svddiv01.gif

 

Notice that although each term is not necessarily linear, the equation as a whole is a linear combination of terms. SVDDIV solves the set of simultaneous equations to yield:

 

image\svddiv02.gif

Remarks:

Consider the problem of fitting a set of data points to a model that is a linear combination of functions, called basis functions:

 

image\svddiv03.gif

 

The basis functions themselves are not required to be linear, only the dependence on the coefficients ck is linear.

 

By employing the method of least squares, we wish to find the coefficients that minimize the sum of the squares of the differences between the data points and the fitted results:

 

image\svddiv04.gif

 

where N is the number of data points and M is the number of basis functions. We define the NxM design matrix:

 

image\svddiv07.gif

 

where the design matrix has the form:

 

image\svddiv08.gif

 

Each column of the design matrix is a basis function evaluated at the known x values.

 

We also define the output vector y of length N as:

 

image\svddiv09.gif

 

With these matrix definitions, we restate the least squares criteria as finding the coefficient vector c that minimizes:

 

image\svddiv24.gif

 

This quantity is minimized to 0 when:

 

image\svddiv13.gif

 

SVDDIV solves the over determined set of simultaneous equations as specified by A and y. Given the matrix equation:

 

A *^ c = y

 

SVDDIV calculates:

 

c = V *^ (1/W) *^ (transpose(U) *^ y)

 

where

 

A = U *^ W *^ transpose(V)

 

as calculated by SVD. By specifying TOL, singular values less than TOL are eliminated and SVDDIV essentially calculates a least squares fit. A typical value for TOL is EPS. See [1] for further details.

 

Although computationally intensive, SVDDIV is particularly useful when the least squares set of equations is ill-conditioned and susceptible to roundoff errors.

 

See LINFIT for a least squares method that uses faster QR decomposition, direct normal equations or SVD and a discussion of the normal equations of the least squares problem.

 

See DADiSP/MatrixXL to significantly optimize SVDDIV.

See Also:

*^ (Matrix Multiply)

\^ (Matrix Solve)

DADiSP/MatrixXL

EPS

LINFIT

MMULT

PINV

POLYFIT

QR

SVD

References:

[1]

Press, Flannery, Teukolsky, Vetterling

 

Numerical Recipes in C

 

Cambridge Press, 1988

 

pp. 407-552