Commit f719b4a3 authored by Henrik Zimmer's avatar Henrik Zimmer

SparseQR now works with Eigen::SparseMatrix and can compute min2norm solutions...

SparseQR now works with Eigen::SparseMatrix and can compute min2norm solutions to underdetermined problem.

git-svn-id: http://www.openflipper.org/svnrepo/CoMISo/trunk@159 1355f012-dd97-4b2f-ae87-10fa9f823a57
parent 77bb4dde
......@@ -336,6 +336,10 @@ endif()
if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_miqp/CMakeLists.txt" )
add_subdirectory (Examples/small_miqp)
endif()
if( EXISTS "${CMAKE_SOURCE_DIR}/Examples/small_sparseqr/CMakeLists.txt" )
add_subdirectory (Examples/small_sparseqr)
endif()
include (ACGCommon)
include (CoMISoExample)
# source code directories
set (directories
.
)
# collect all header and source files
acg_append_files (headers "*.hh" ${directories})
acg_append_files (sources "*.cc" ${directories})
# remove template cc files from source file list
acg_drop_templates (sources)
if (WIN32)
acg_add_executable (small_sparseqr WIN32 ${sources} ${headers} )
elseif (APPLE)
# generate bundle on mac
acg_add_executable (small_sparseqr MACOSX_BUNDLE ${sources} ${headers} )
else ()
acg_add_executable (small_sparseqr ${sources} ${headers} )
endif ()
# enable rpath linking
set_target_properties(small_sparseqr PROPERTIES INSTALL_RPATH_USE_LINK_PATH 1)
target_link_libraries (small_sparseqr
CoMISo
${COMISO_LINK_LIBRARIES}
)
if (APPLE)
# create bundle in "Build" directory and set icon
# no install needed here, because the whole bundle will be installed in the
# toplevel CMakeLists.txt
set_target_properties (
small_sparseqr PROPERTIES
RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/Build"
MACOSX_BUNDLE_INFO_STRING "CoMISo small_sparseqr"
)
endif ()
/*===========================================================================*\
* *
* CoMISo *
* Copyright (C) 2008-2009 by Computer Graphics Group, RWTH Aachen *
* www.rwth-graphics.de *
* *
*---------------------------------------------------------------------------*
* This file is part of CoMISo. *
* *
* CoMISo 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. *
* *
* CoMISo 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 CoMISo. If not, see <http://www.gnu.org/licenses/>. *
* *
\*===========================================================================*/
#include <CoMISo/Config/config.hh>
#include <CoMISo/Utils/StopWatch.hh>
#include <vector>
#include <cstdlib>
#include <Eigen/Sparse>
#include <CoMISo/Solver/SparseQRSolver.hh>
#include <CoMISo/Solver/Eigen_Tools.hh>
//------------------------------------------------------------------------------------------------------
// Example main
int main(void)
{
std::cout << "---------- 0) Using Sparse QR for solving underdetermined equations and computing Null spaces " << std::endl;
typedef Eigen::SparseMatrix< double > SpMatrix;
typedef Eigen::MatrixXd DenMatrix;
typedef Eigen::Triplet< double > Triplet;
int dimr(4+1);
int dimc(4+2);
std::cout << "---------- 1) Creating matrix " << std::endl;
std::vector< Triplet > triplets;
for( int i = 0; i < dimc*dimr/2; ++i)
{
int x( rand()%(dimr-1));
int y( rand()%dimc);
double val( rand()%10);
//std::cerr << " setting (" << x << ", " << y << ") to " << val << std::endl;
triplets.push_back( Triplet( x, y, val));
}
SpMatrix A(dimr,dimc);
A.setFromTriplets(triplets.begin(), triplets.end());
std::cerr << DenMatrix(A) << std::endl;
int m = dimr;
int n = dimc;
if( m < n )
{
std::swap( m,n);
std::cerr << " ... m < n -> form transposed ..." << std::endl;
A = SpMatrix(A.transpose());
// test make also row -rank-deficinet
A.middleCols(n-1,1) = A.middleCols(0,1);
A.middleCols(0,1) = A.middleCols(n-2,1);
std::cerr << DenMatrix(A) << std::endl;
}
std::cerr << " ... m = " << m << "; n = " << n << std::endl;
std::cerr << std::endl;
std::cout << "---------- 2) Sparse QR " << std::endl;
COMISO::SparseQRSolver spqr;
SpMatrix Q,R;
std::vector< size_t > P;
int rank = spqr.factorize_system_eigen( A, Q, R, P);
int nullity(dimc-rank);
// setup permutation matrix
SpMatrix Pm( n, n);
if( P.size() != 0)
{
for( size_t i = 0; i < P.size(); ++i)
{
Pm.coeffRef( i, P[i]) = 1;
}
}
std::cout << "---------- 3) Result " << std::endl;
std::cerr << " Q " << std::endl << DenMatrix(Q) << std::endl;
std::cerr << " R " << std::endl << DenMatrix(R) << std::endl;
std::cerr << " P " << std::endl << P << std::endl;
std::cerr << " P matrix " << std::endl << DenMatrix(Pm) << std::endl;
std::cerr << " Rank " << rank << std::endl;
std::cerr << " Nullity " << nullity << std::endl;
// extract nullspace
SpMatrix NullSpace( Q.middleCols( std::max( 0, m-nullity), nullity));
std::cerr << " Nullspace " << std::endl << DenMatrix(NullSpace) << std::endl;
// non nullspace part of R
//// assuming superflous column in R is the last (if A is also row deficient)
//SpMatrix Rtmp(R.middleCols(0,std::min(n,n-(n-rank))).transpose());
//SpMatrix R1( R.transpose().middleCols(0, m-nullity));
SpMatrix Rtmp(R.transpose());
SpMatrix R1t( Rtmp.middleCols(0,m-nullity));
SpMatrix R1( R1t.transpose());
std::cerr << " Non-Nullspace R " << std::endl << DenMatrix(R1) << std::endl;
std::cout << "---------- 4) Verification " << std::endl;
SpMatrix reconstructedA(Q*R*Pm.transpose());
std::cerr << " Q orthogonal? \t " << ((fabs((Q.transpose()*Q).squaredNorm()-m) < 1e-8)?"yes":"no") << std::endl;
std::cerr << " A = QR? \t " << (((reconstructedA-A).squaredNorm() < 1e-8)? "yes":"no") << std::endl;
std::cerr << std::endl << std::endl;
std::cout << "---------- 5) Solving Ax=b (with x without nullspace component)" << std::endl;
// NOTE: A was transposed above to be m>n
SpMatrix b(n,1);
SpMatrix x(m,1);
for( int i = 0; i < n; ++i)
b.coeffRef(i,0) = rand()%10;
std::cerr << " ... System Ax = b .. \n";
std::cerr << " A " << std::endl << DenMatrix(A.transpose()) << " x " << std::endl << DenMatrix(x) << " b " << std::endl << DenMatrix(b) << std::endl;
std::cout << "---------- 5.1) test: solve using sparse QR solving .." << std::endl;
SpMatrix At(A.transpose());
spqr.solve_system_eigen( At, b, x);
std::cerr << " ... solution x .. " << std::endl;
std::cerr << DenMatrix(x) << std::endl;
std::cerr << " ... test: is a solution ? " << (((A.transpose()*x-b).squaredNorm()<1e-8)?"yes":"no") << std::endl;
std::cerr << " ... test: has nullspace component ? " << ((x.transpose()*NullSpace).squaredNorm()<1e-8?"yes":"no") << std::endl;
std::cerr << " ... Nullspace projections : " << (x.transpose()*NullSpace) << std::endl;
std::cout << "---------- 5.2) test: solve without nullspace .." << std::endl;
SpMatrix Atnull(At);
SpMatrix bnull(b);
SpMatrix xnull(m,1);
spqr.solve_system_eigen_min2norm( Atnull, bnull, xnull);
std::cerr << " ... solution x .. " << std::endl;
std::cerr << DenMatrix(xnull) << std::endl;
std::cerr << " ... test: is a solution ? " << (((A.transpose()*xnull-bnull).squaredNorm()<1e-8)?"yes":"no") << std::endl;
std::cerr << " ... test: has nullspace component ? " << ((xnull.transpose()*NullSpace).squaredNorm()<1e-8?"yes":"no") << std::endl;
std::cerr << " ... Nullspace projections : " << (xnull.transpose()*NullSpace) << std::endl;
return 0;
}
......@@ -276,6 +276,367 @@ bool is_symmetric( const MatrixT& _A)
return symmetric;
}
//-----------------------------------------------------------------------------
template< class Eigen_MatrixT, class IntT >
void permute( const Eigen_MatrixT& _QR, const std::vector< IntT>& _Pvec, Eigen_MatrixT& _A)
{
typedef typename Eigen_MatrixT::Scalar Scalar;
int m = _QR.innerSize();
int n = _QR.outerSize();
if( _Pvec.size() == 0)
{
_A = _QR;
return;
}
if( _Pvec.size() != (size_t)_QR.cols() && _Pvec.size() != 0)
{
std::cerr << __FUNCTION__ << " wrong size of permutation vector, should have #cols length (or zero)" << std::endl;
}
// build sparse permutation matrix
typedef Eigen::Triplet< Scalar > Triplet;
std::vector< Triplet > triplets;
triplets.reserve(_QR.nonZeros());
_A = Eigen_MatrixT( m, n);
typedef typename Eigen_MatrixT::InnerIterator It;
for( int c = 0; c < n; ++c) // cols
{
for( It it(_QR,c); it; ++it) // rows
{
int r(it.index());
Scalar val(it.value());
int newcol( _Pvec[c]);
triplets.push_back( Triplet( r, newcol, val));
}
}
_A.setFromTriplets( triplets.begin(), triplets.end());
}
//-----------------------------------------------------------------------------
#ifndef COMISO_NCHOLMOD
/// Eigen to Cholmod_sparse interface
template<class MatrixT>
void cholmod_to_eigen( const cholmod_sparse& _AC, MatrixT& _A)
{
// initialize dimensions
typedef typename MatrixT::Scalar Scalar;
typedef Eigen::Triplet< Scalar > Triplet;
size_t nzmax( _AC.nzmax);
std::cerr << __FUNCTION__ << " row " << _AC.nrow << " col " << _AC.ncol << " stype " << _AC.stype << std::endl;
_A = MatrixT(_AC.nrow, _AC.ncol);
std::vector< Triplet > triplets;
triplets.reserve(nzmax);
if(!_AC.packed)
{
std::cerr << "Warning: " << __FUNCTION__ << " does not support unpacked matrices yet!!!" << std::endl;
return;
}
// Pointer to data
double* X((double*)_AC.x);
// complete matrix stored
if(_AC.stype == 0)
{
// which integer type?
if(_AC.itype == CHOLMOD_LONG)
{
UF_long* P((UF_long*)_AC.p);
UF_long* I((UF_long*)_AC.i);
for(UF_long i=0; i<(UF_long)_AC.ncol; ++i)
for(UF_long j= P[i]; j< P[i+1]; ++j)
//_A( I[j], i) += X[j]; // += really needed?
triplets.push_back( Triplet( I[j], i, X[j]));
}
else
{
int* P((int*)_AC.p);
int* I((int*)_AC.i);
for(int i=0; i<(int)_AC.ncol; ++i)
for(int j= P[i]; j< P[i+1]; ++j)
triplets.push_back( Triplet( I[j], i, X[j]));
//_A( I[j], i) += X[j];
}
}
else // only upper or lower diagonal stored
{
// which integer type?
if(_AC.itype == CHOLMOD_LONG)
{
UF_long* P((UF_long*)_AC.p);
UF_long* I((UF_long*)_AC.i);
for(UF_long i=0; i<(UF_long)_AC.ncol; ++i)
for(UF_long j=P[i]; j<P[i+1]; ++j)
{
//_A(I[j], i) += X[j];
triplets.push_back( Triplet( I[j], i, X[j]));
// add up symmetric part
if( I[j] != i)
triplets.push_back( Triplet( i, I[j], X[j]));
//_A(i,I[j]) += X[j];
}
}
else
{
int* P((int*)_AC.p);
int* I((int*)_AC.i);
for(int i=0; i<(int)_AC.ncol; ++i)
for(int j=P[i]; j<P[i+1]; ++j)
{
//_A(I[j], i) += X[j];
triplets.push_back( Triplet( I[j], i, X[j]));
// add up symmetric part
if( I[j] != i)
// _A(i,I[j]) += X[j];
triplets.push_back( Triplet( i, I[j], X[j]));
}
}
}
_A.setFromTriplets( triplets.begin(), triplets.end());
}
/// GMM to Cholmod_sparse interface
template<class MatrixT>
void eigen_to_cholmod( const MatrixT& _A, cholmod_sparse* &_AC, cholmod_common* _common, int _sparsity_type, bool _long_int)
{
/* _sparsity_type
* 0: matrix is "unsymmetric": use both upper and lower triangular parts
* (the matrix may actually be symmetric in pattern and value, but
* both parts are explicitly stored and used). May be square or
* rectangular.
* >0: matrix is square and symmetric, use upper triangular part.
* Entries in the lower triangular part are ignored.
* <0: matrix is square and symmetric, use lower triangular part.
* Entries in the upper triangular part are ignored. */
int m = _A.innerSize();
int n = _A.outerSize();
// get upper or lower
char uplo = 'c';
if(_sparsity_type < 0) uplo = 'l';
if(_sparsity_type > 0) uplo = 'u';
if( _long_int) // long int version
{
std::vector<double> values;
std::vector<UF_long> rowind;
std::vector<UF_long> colptr;
// get data of gmm matrix
COMISO_EIGEN::get_ccs_symmetric_data( _A, uplo, values, rowind, colptr);
// allocate cholmod matrix
_AC = cholmod_l_allocate_sparse(m,n,values.size(),true,true,_sparsity_type,CHOLMOD_REAL, _common);
// copy data to cholmod matrix
for(UF_long i=0; i<(UF_long)values.size(); ++i)
{
((double*) (_AC->x))[i] = values[i];
((UF_long*)(_AC->i))[i] = rowind[i];
}
for(UF_long i=0; i<(UF_long)colptr.size(); ++i)
((UF_long*)(_AC->p))[i] = colptr[i];
}
else // int version
{
std::vector<double> values;
std::vector<int> rowind;
std::vector<int> colptr;
// get data of gmm matrix
COMISO_EIGEN::get_ccs_symmetric_data( _A, uplo, values, rowind, colptr);
// allocate cholmod matrix
_AC = cholmod_allocate_sparse(m,n,values.size(),true,true,_sparsity_type,CHOLMOD_REAL, _common);
// copy data to cholmod matrix
for(unsigned int i=0; i<values.size(); ++i)
{
((double*)(_AC->x))[i] = values[i];
((int*) (_AC->i))[i] = rowind[i];
}
for(unsigned int i=0; i<colptr.size(); ++i)
((int*)(_AC->p))[i] = colptr[i];
}
}
/*
/// Eigen to Cholmod_dense interface
template<class MatrixT>
void cholmod_to_eigen_dense( const cholmod_dense& _AC, MatrixT& _A)
{
// initialize dimensions
typedef typename MatrixT::Scalar Scalar;
typedef Eigen::Triplet< Scalar > Triplet;
size_t nzmax( _AC.nzmax);
_A = MatrixT(_AC.nrow, _AC.ncol);
std::vector< Triplet > triplets;
triplets.reserve(nzmax);
if(!_AC.packed)
{
std::cerr << "Warning: " << __FUNCTION__ << " does not support unpacked matrices yet!!!" << std::endl;
return;
}
// Pointer to data
double* X((double*)_AC.x);
// complete matrix stored
if(_AC.stype == 0)
{
// which integer type?
if(_AC.itype == CHOLMOD_LONG)
{
UF_long* P((UF_long*)_AC.p);
UF_long* I((UF_long*)_AC.i);
for(UF_long i=0; i<(UF_long)_AC.ncol; ++i)
for(UF_long j= P[i]; j< P[i+1]; ++j)
//_A( I[j], i) += X[j]; // += really needed?
triplets.push_back( Triplet( I[j], i, X[j]));
}
else
{
int* P((int*)_AC.p);
int* I((int*)_AC.i);
for(int i=0; i<(int)_AC.ncol; ++i)
for(int j= P[i]; j< P[i+1]; ++j)
triplets.push_back( Triplet( I[j], i, X[j]));
//_A( I[j], i) += X[j];
}
}
else // only upper or lower diagonal stored
{
// which integer type?
if(_AC.itype == CHOLMOD_LONG)
{
UF_long* P((UF_long*)_AC.p);
UF_long* I((UF_long*)_AC.i);
for(UF_long i=0; i<(UF_long)_AC.ncol; ++i)
for(UF_long j=P[i]; j<P[i+1]; ++j)
{
//_A(I[j], i) += X[j];
triplets.push_back( Triplet( I[j], i, X[j]));
// add up symmetric part
if( I[j] != i)
triplets.push_back( Triplet( i, I[j], X[j]));
//_A(i,I[j]) += X[j];
}
}
else
{
int* P((int*)_AC.p);
int* I((int*)_AC.i);
for(int i=0; i<(int)_AC.ncol; ++i)
for(int j=P[i]; j<P[i+1]; ++j)
{
//_A(I[j], i) += X[j];
triplets.push_back( Triplet( I[j], i, X[j]));
// add up symmetric part
if( I[j] != i)
// _A(i,I[j]) += X[j];
triplets.push_back( Triplet( i, I[j], X[j]));
}
}
}
_A.setFromTriplets( triplets.begin(), triplets.end());
}
/// GMM to Cholmod_sparse interface
template<class MatrixT>
void eigen_to_cholmod_dense( const MatrixT& _A, cholmod_dense* &_AC, cholmod_common* _common, bool _long_int)
{
int m = _A.innerSize();
int n = _A.outerSize();
// allocate cholmod matrix
_AC = cholmod_l_allocate_sparse(m,n,values.size(),true,true,_sparsity_type,CHOLMOD_REAL, _common);
_AC = cholmod_l_allocate_dense (m,n,n, xtype, cc)
if( _long_int) // long int version
{
std::vector<double> values;
std::vector<UF_long> rowind;
std::vector<UF_long> colptr;
// get data of gmm matrix
COMISO_EIGEN::get_ccs_symmetric_data( _A, uplo, values, rowind, colptr);
// allocate cholmod matrix
_AC = cholmod_l_allocate_sparse(m,n,values.size(),true,true,_sparsity_type,CHOLMOD_REAL, _common);
// copy data to cholmod matrix
for(UF_long i=0; i<(UF_long)values.size(); ++i)
{
((double*) (_AC->x))[i] = values[i];
((UF_long*)(_AC->i))[i] = rowind[i];
}
for(UF_long i=0; i<(UF_long)colptr.size(); ++i)
((UF_long*)(_AC->p))[i] = colptr[i];
}
else // int version
{
std::vector<double> values;
std::vector<int> rowind;
std::vector<int> colptr;
// get data of gmm matrix
COMISO_EIGEN::get_ccs_symmetric_data( _A, uplo, values, rowind, colptr);
// allocate cholmod matrix
_AC = cholmod_allocate_sparse(m,n,values.size(),true,true,_sparsity_type,CHOLMOD_REAL, _common);
// copy data to cholmod matrix
for(unsigned int i=0; i<values.size(); ++i)
{
((double*)(_AC->x))[i] = values[i];
((int*) (_AC->i))[i] = rowind[i];
}
for(unsigned int i=0; i<colptr.size(); ++i)
((int*)(_AC->p))[i] = colptr[i];
}
}*/
#endif
//=============================================================================
} // namespace COMISO
......
......@@ -34,6 +34,8 @@
#include <algorithm>
#include <limits>
#include <cmath>
#include <Eigen/Eigen>
#include <Eigen/Sparse>
//#ifdef COMISO_Eigen3_AVAILABLE
//#include <Eigen/Eigen>
......@@ -87,6 +89,23 @@ void inspect_matrix( const MatrixT& _A);
template<class MatrixT>
bool is_symmetric( const MatrixT& _A);
template< class Eigen_MatrixT, class IntT >
void permute( const Eigen_MatrixT& _QR, const std::vector< IntT>& _Pvec, Eigen_MatrixT& _A);
#ifndef COMISO_NCHOLMOD
/// Eigen to Cholmod_sparse interface
template<class MatrixT>
void cholmod_to_eigen( const cholmod_sparse& _AC, MatrixT& _A);