/*****************************************************************************
*                                                                            *
*   SVDDIV.SPL    Copyright (C) 1999-2018 DSP Development Corporation        *
*                               All Rights Reserved                          *
*                                                                            *
*   Author:       Randy Race                                                 *
*                                                                            *
*   Synopsis:     Solves A *^ c = x using singular value decomposition       *
*                                                                            *
*   Revisions:    18 Jun 1999  RRR  Creation                                 *
*                 21 Mar 2000  RRR  use input matrix column size             *
*                 18 Jan 2010  RRR  use faster economy SVD mode              *
*                 26 Sep 2018  RRR  "econ" fix                               *
*                                                                            *
*****************************************************************************/

#include 

#if @HELP_SVDDIV

    SVDDIV

    Purpose: Solves for c in A *^ c = y using singular value decomposition

    Syntax:  SVDDIV(A, y, tol)

             A   - input array, the basis function design matrix

             y   - input series, the Y output series

             tol - optional real, threshold at which to zero out singular
                   values. If not specified, all singular values are used.

    Returns: A series or array


    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:

             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: svddiv(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. 
             SVDDIV 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. 
             This quantity is minimized to 0 when:

                          A *^ c = y

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

                          A *^ c = y

             SVDDIV calculates:

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

              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.

             See LINFIT for a least squares method that uses faster QR
             decomposition, direct normal equations or SVD.

    See Also:
             *^
             LINFIT
             QR
             SVD

    References:

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



/* solve A *^ c = y using singular value decomposition */
svddiv(A, y, tol)
{
        local m, n, c, u, w, v, sing_edit = TRUE, coefstd;

        if (argc < 3)
        {
                if (argc < 2) error("svddiv - two input matrices required");
                sing_edit = FALSE;
        }

        /* get num rows & columns */
        (m, n) = size(A);

        /* SVD matrices - use economy sized matrix where u is NxM */
        (u, w, v) = svd(A, "econ");

        /* get diagonal elements */
        w = diag(extract(w, 1, n));

        /* convert to diagonal matrix, edit singular values if necessary */
        if (sing_edit)
        {
                /* edit singular values, if w < tol then 1/w = 0 */
                w = diag((1 / w) * (w > tol));
        }
        else   
        {
                /* use inverse of computed singular values */
                w = diag(1 / w);
        }

        /* solve the system */
        c = v *^ (w *^ (transposeconj(u) *^ y));

        if (outargc > 1)
        {
                /* return coefficient stdev, note that w already inversed */
                coefstd = sqrt(transpose(colsum((transpose(v) * diag(w)) ^ 2)));

                return(c, coefstd);
        }
        else
        {
                return(c);
        }
}