/home/jgoppert/Projects/uvsim/uvsim/src/test/pinv.cc

Go to the documentation of this file.

/*
 * pinv.cc
 * Copyright (C) James Goppert 2009 <jgoppert@purdue.edu>
 *
 * pinv.cc is free software: you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * pniv.cc 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 GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License along
 * with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <iostream>
#include <sstream>
#include <iomanip>
#include <term.h>
#include <unistd.h>
#include <cmath>
#include "math.h"
#include <limits>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
#include <boost/numeric/bindings/traits/ublas_vector.hpp>
#include <boost/numeric/bindings/lapack/gesvd.hpp>
#include <boost/numeric/bindings/lapack/workspace.hpp>
#include <string>
#include <algorithm>

using namespace boost::numeric::ublas;
using namespace boost::numeric::bindings;

double max(vector<double> v);
matrix<double> pinv(matrix<double> A);

int main() 
{
        matrix<double> A(3,3);
        A(0,0) = 1;
        A(1,2) = 2;
        A(2,2) = .2;
        A(2,1) = .5;
        std::cout << "A: " << A << std::endl;
        std::cout << "pinv(A): " << pinv(A) << std::endl;
}

double max(vector<double> v)
{
    double max = v(0);
    for (int i=1;i<v.size();i++)
    {
        if (v(i) > max) max = v(i);
    }
    return max;
}

matrix<double> pinv(matrix<double> A)
{
    matrix<double, column_major> Ac = A;
    int m = A.size1(), n = A.size2();
    int maxDim = std::max(m,n), minDim = std::min(m,n);
    vector<double> Sig(minDim);
    matrix<double, column_major> U(m,m), SigI(n,m), VT(n,n);
    std::cout << "Performing SVD" << std::endl;
    lapack::gesvd('A','A', Ac, Sig, U, VT);

    // Computer pseudo inverse of Sigma
    double tol = std::numeric_limits<double>::epsilon()*maxDim*max(Sig);
    //std::cout << "epsilon: " << std::numeric_limits<double>::epsilon() << std::endl;
    //std::cout << "maxDim: " << maxDim << std::endl;
    //std::cout << "maxSig: " << max(Sig) << std::endl;
    //std::cout << "tol: " << tol << std::endl;

    SigI = zero_matrix<double>(n,m);
    for (int i=0;i<minDim;i++)
    {
        if (Sig(i) > tol)
        {
            SigI(i,i) = 1/Sig(i);
        }
        else
        {
            SigI(i,i) = 0;
        }
    }

    // Debug
    //std::cout << "U: " << U << std::endl;
    //std::cout << "Sig: " << Sig << std::endl;
    //std::cout << "SigI: " << SigI << std::endl;
    //std::cout << "V: " << trans(VT) << std::endl;

    // Computer pseudo inverse of A
    matrix<double> temp = prod(trans(VT),SigI);
    return prod(temp,trans(U));
}

// vim:ts=4:sw=4

---