/*****************************************************************************
*                                                                            *
*   LINFIT.SPL    Copyright (C) 2010 DSP Development Corporation             *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:       Randy Race                                                 *
*                                                                            *
*   Synopsis:     Solves A *^ y = c using QR, SVD or normal equations        *
*                                                                            *
*   Revisions:    11 Jan 2010  RRR  Creation                                 *
*                                                                            *
*****************************************************************************/

#include 

#if @HELP_LINFIT

    LINFIT

    Purpose: Performs general least squares fitting of basis functions 
             using QR decomposition, SVD or normal equations.

    Syntax:  LINFIT(A, y, method)

                  A - input array, the basis function design matrix

                  y - input series, the Y output series

             method - a string, the fitting method

                        "qr"   - QR decomposition (default)
                        "svd"  - singular value decomposition
                        "norm" - normal equations

    Returns: A series, the basis function coefficients.

    Example:
             x = {-1, 1, 1.5, 3};
             y = {6, 2, 2.25, 6};

             c1 = polyfit(y, x, 2);
             c2 = linfit(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, LINFIT produces
             the same result. The fitted equation becomes:

             y(x) = 3 -2 * x + x^2

    Example:
             W1: gnorm(100,.01)
             W2: 5*w1 + 3*sin(w1) - 2*w1^3
             W3: ravel(w1, sin(w1), w1^3)
             W4: linfit(W3, w2)

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

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

             y(x) = c1 * x + c2 * sin(x) + c3 * x^3 

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

             c1 ==  5
             c2 ==  3
             c3 == -2
 
    Remarks:

             For the least squares problem, we wish to find the coefficients
             c that minimize:

                          |A *^c - y|^2

             Where A is the design matrix and y is the output vector. 

             By taking the partial derivatives with respect to each
             coefficient, the least squares problem becomes:

                         A' *^ A = A' *^ y

             The equations or known as the normal equations for the least
             squares problem.

             LINFIT(A, y, "norm") solves for the normal equations
             directly.  However, the normal equations can often be
             ill-conditioned and solving for the coefficients directly
             can be susceptible to roundoff errors.

             To investigate an alternative solution to the least
             squares problem, notice that the least squares merit function
             is minimized to 0 when:

                          A *^ c = y

             LINFIT(A, y) or LINFIT(A, y, "qr") solves the set of
             simultaneous equations by QR decomposition where A is
             an N x M array, y is a vector of length N, N is the number
             of data points and M is the number of basis functions.

             Using QR decomposition, matrix A is expressed as the
             product of an N x N orthogonal matrix Q and a N x M matrix
             R where R is an upper triangular matrix such that the last
             N-M rows contain zeros.  

             LINFIT calculates:

                Q' *^ Q *^ R *^ c = Q' *^ y

                R *^ c = Q' *^ y

            Since R is an upper triangular matrix, the system can be solved
            with simple back substitution.

            Singular value decomposition decomposes the design matrix A
            into three matrix components, U, W and V such that:

                A = U *^ W *^ V'

            where U is an orthonormal matrix that contains the left
            singular values, W is a diagonal matrix that contains the
            weighting factor or singular values and V is a orthonormal
            matrix that contains the right singular values.

            With this decomposition, the least squares problem can be
            solved with:

                c = V *^ diag(diag(W)) *^ U' *^ y

            LINFIT(A, y, "svd") uses singular value decomposition to
            solve the over-determined set of simultaneous equations as
            specified by A and y.  The method is usually slower than QR
            decomposition, but can be more stable for severely
            ill-conditioned systems.

            See [1] for further details.

    See Also:
             \^
             POLYFIT
             QR
             SVD
             SVDDIV

    References:

             [1] Press, Flannery, Teukolsky, Vetterling
                 Numerical Recipes in C
                 2nd Edition
                 Cambridge Press, 1988, 1992
                 pp. 676-680
#endif



/* solve A *^ c = y using QR decomposition */
linfit(A, y, method)
{
        local c;

        if (argc < 3)
        {
                if (argc < 2) error("linfit - two input matrices required");

                method = "qr";
        }

        switch (tolower(method))
        {
                case "qr":
                default:
                        /* solve overdetermined system via QR */
                        c = A \^ y;
                        break;

                case "svd":
                        /* solve overdetermined system via SVD */
                        c = svddiv(A, y);
                        break;

                case "norm":
                        /* solve via normal equations */
                        c = crossprod(A) \^ crossprod(A, y);
                        break;
        }

        return(c);
}