svd.cpp

有很多的函数库

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2003 Neil Firth

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <[e-mail]>. The license is also available online at
 <http://quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.

        Adapted from the TNT project
        http://math.nist.gov/tnt/download.html

        This software was developed at the National Institute of Standards
        and Technology (NIST) by employees of the Federal Government in the
        course of their official duties. Pursuant to title 17 Section 105
        of the United States Code this software is not subject to copyright
        protection and is in the public domain. NIST assumes no responsibility
        whatsoever for its use by other parties, and makes no guarantees,
        expressed or implied, about its quality, reliability, or any other
        characteristic.

        We would appreciate acknowledgement if the software is incorporated in
        redistributable libraries or applications.
*/


#include <ql/math/matrixutilities/svd.hpp>

namespace QuantLib {

    namespace {

        /*  returns hypotenuse of real (non-complex) scalars a and b by
            avoiding underflow/overflow
            using (a * sqrt( 1 + (b/a) * (b/a))), rather than
            sqrt(a*a + b*b).
        */
        Real hypot(const Real &a, const Real &b) {
            if (a == 0) {
                return std::fabs(b);
            } else {
                Real c = b/a;
                return std::fabs(a) * std::sqrt(1 + c*c);
            }
        }

    }


    SVD::SVD(const Matrix& M) {

        using std::swap;

        Matrix A;

        /* The implementation requires that rows > columns.
           If this is not the case, we decompose M^T instead.
           Swapping the resulting U and V gives the desired
           result for M as

           M^T = U S V^T           (decomposition of M^T)

           M = (U S V^T)^T         (transpose)

           M = (V^T^T S^T U^T)     ((AB)^T = B^T A^T)

           M = V S U^T             (idempotence of transposition,
                                    symmetry of diagonal matrix S)

        */

        if (M.rows() >= M.columns()) {
            A = M;
            transpose_ = false;
        } else {
            A = transpose(M);
            transpose_ = true;
        }

        m_ = A.rows();
        n_ = A.columns();

        // we're sure that m_ >= n_

        s_ = Array(n_);
        U_ = Matrix(m_,n_, 0.0);
        V_ = Matrix(n_,n_);
        Array e(n_);
        Array work(m_);
        Integer i, j, k;

        // Reduce A to bidiagonal form, storing the diagonal elements
        // in s and the super-diagonal elements in e.

        Integer nct = std::min(m_-1,n_);
        Integer nrt = std::max(0,n_-2);
        for (k = 0; k < std::max(nct,nrt); k++) {
            if (k < nct) {

                // Compute the transformation for the k-th column and
                // place the k-th diagonal in s[k].
                // Compute 2-norm of k-th column without under/overflow.
                s_[k] = 0;
                for (i = k; i < m_; i++) {
                    s_[k] = hypot(s_[k],A[i][k]);
                }
                if (s_[k] != 0.0) {
                    if (A[k][k] < 0.0) {
                        s_[k] = -s_[k];
                    }
                    for (i = k; i < m_; i++) {
                        A[i][k] /= s_[k];
                    }
                    A[k][k] += 1.0;
                }
                s_[k] = -s_[k];
            }
            for (j = k+1; j < n_; j++) {
                if ((k < nct) && (s_[k] != 0.0))  {

                    // Apply the transformation.

                    Real t = 0;
                    for (i = k; i < m_; i++) {
                        t += A[i][k]*A[i][j];
                    }
                    t = -t/A[k][k];
                    for (i = k; i < m_; i++) {
                        A[i][j] += t*A[i][k];
                    }
                }

                // Place the k-th row of A into e for the
                // subsequent calculation of the row transformation.

                e[j] = A[k][j];
            }
            if (k < nct) {

                // Place the transformation in U for subsequent back
                // multiplication.

                for (i = k; i < m_; i++) {
                    U_[i][k] = A[i][k];
                }
            }
            if (k < nrt) {

                // Compute the k-th row transformation and place the
                // k-th super-diagonal in e[k].
                // Compute 2-norm without under/overflow.
                e[k] = 0;
                for (i = k+1; i < n_; i++) {
                    e[k] = hypot(e[k],e[i]);
                }
                if (e[k] != 0.0) {
                    if (e[k+1] < 0.0) {
                        e[k] = -e[k];
                    }
                    for (i = k+1; i < n_; i++) {
                        e[i] /= e[k];
                    }
                    e[k+1] += 1.0;
                }
                e[k] = -e[k];
                if ((k+1 < m_) & (e[k] != 0.0)) {

                    // Apply the transformation.

                    for (i = k+1; i < m_; i++) {
                        work[i] = 0.0;
                    }
                    for (j = k+1; j < n_; j++) {
                        for (i = k+1; i < m_; i++) {
                            work[i] += e[j]*A[i][j];
                        }
                    }
                    for (j = k+1; j < n_; j++) {
                        Real t = -e[j]/e[k+1];
                        for (i = k+1; i < m_; i++) {
                            A[i][j] += t*work[i];
                        }
                    }
                }

                // Place the transformation in V for subsequent
                // back multiplication.

                for (i = k+1; i < n_; i++) {
                    V_[i][k] = e[i];
                }
            }
        }

        // Set up the final bidiagonal matrix or order n.

        if (nct < n_) {
            s_[nct] = A[nct][nct];
        }
        if (nrt+1 < n_) {
            e[nrt] = A[nrt][n_-1];
        }
        e[n_-1] = 0.0;

        // generate U

        for (j = nct; j < n_; j++) {
            for (i = 0; i < m_; i++) {
                U_[i][j] = 0.0;
            }
            U_[j][j] = 1.0;
        }
        for (k = nct-1; k >= 0; --k) {
            if (s_[k] != 0.0) {
                for (j = k+1; j < n_; ++j) {
                    Real t = 0;
                    for (i = k; i < m_; i++) {
                        t += U_[i][k]*U_[i][j];
                    }
                    t = -t/U_[k][k];
                    for (i = k; i < m_; i++) {
                        U_[i][j] += t*U_[i][k];
                    }
                }
                for (i = k; i < m_; i++ ) {
                    U_[i][k] = -U_[i][k];
                }
                U_[k][k] = 1.0 + U_[k][k];
                for (i = 0; i < k-1; i++) {
                    U_[i][k] = 0.0;
                }
            } else {
                for (i = 0; i < m_; i++) {
                    U_[i][k] = 0.0;
                }
                U_[k][k] = 1.0;
            }
        }

        // generate V

        for (k = n_-1; k >= 0; --k) {
            if ((k < nrt) & (e[k] != 0.0)) {
                for (j = k+1; j < n_; ++j) {
                    Real t = 0;
                    for (i = k+1; i < n_; i++) {
                        t += V_[i][k]*V_[i][j];
                    }
                    t = -t/V_[k+1][k];
                    for (i = k+1; i < n_; i++) {
                        V_[i][j] += t*V_[i][k];
                    }
                }
            }
            for (i = 0; i < n_; i++) {
                V_[i][k] = 0.0;
            }
            V_[k][k] = 1.0;
        }

        // Main iteration loop for the singular values.