9

Does boost have one? Where A, y and x is a matrix (sparse and can be very large) and vectors respectively. Either y or x can be unknown.

I can't seem to find it here: http://www.boost.org/doc/libs/1_39_0/libs/numeric/ublas/doc/index.htm

Eric Wilson
  • 57,719
  • 77
  • 200
  • 270
neversaint
  • 60,904
  • 137
  • 310
  • 477
  • 2
    Since I know nothing about Boost or C++, I'll just comment that there is a solution if and only if A is an invertible matrix. If A is not invertible, the least squares solution x = (A^T A)^(-1) A^T y is the closest thing. – Eric Wilson Aug 04 '09 at 02:09
  • 4
    that's generally not the way matrix equations are solved (inverses are generally bad, both for numerical stability and for speed), rather you'll see QR or LU factorization followed by back-substitution. – Jason S Aug 04 '09 at 13:38

6 Answers6

19

yes, you can solve linear equations with boost's ublas library. Here is one short way using LU-factorize and back-substituting to get the inverse:

using namespace boost::ublas;

Ainv = identity_matrix<float>(A.size1());
permutation_matrix<size_t> pm(A.size1());
lu_factorize(A,pm)
lu_substitute(A, pm, Ainv);

So to solve a linear system Ax=y, you would solve the equation trans(A)Ax=trans(A)y by taking the inverse of (trans(A)A)^-1 to get x: x = (trans(A)A)^-1Ay.

Inverse
  • 4,408
  • 2
  • 26
  • 35
  • 12
    If all you need is a solution for Ax=y, just use `permutation_matrix pm(A.size1()); lu_factorize(A, pm); lu_substitute(A, pm, y);` and `y` is replaced with the solution. – Joey Jul 10 '14 at 17:24
  • Shouldn't this be the equation for x? `x = (trans(A) * A)^-1 * trans(A) * y` – Synck Jul 19 '18 at 11:14
5

Linear solvers are generally part of the LAPACK library which is a higher level extension of the BLAS library. If you are on Linux, the Intel MKL has some good solvers, optimized both for dense and sparse matrices. If you are on windows, MKL has a one month trial for free... and to be honest I haven't tried any of the other ones out there. I know the Atlas package has a free LAPACK implementation but not sure how hard it is to get running on windows.

Anyways, search around for a LAPACK library which works on your system.

DeusAduro
  • 5,971
  • 5
  • 29
  • 36
  • 4
    Note: LAPACK is not a sparse solver, so it does not store sparse matrices very efficiently, nor does it solve sparse systems particularly efficiently. The Intel MKL does contain sparse solvers ( http://software.intel.com/sites/products/collateral/hpc/mkl/mkl_brief.pdf ), in particular the PARDISO direct sparse solver ( http://www.pardiso-project.org ). A good overview of dense and sparse numerical linear algebra software can be found at http://www.netlib.org/utk/people/JackDongarra/la-sw.html – las3rjock Aug 08 '09 at 12:16
  • Indeed it does, and unless you are building a commercial product I would suggest forgoing the MKL library and just getting PARDISO, you'll save money and save having to deal with licensing issues. – DeusAduro Aug 08 '09 at 19:15
4

One of the best solvers for Ax = b, when A is sparse, is Tim Davis's UMFPACK

UMFPACK computes a sparse LU decomposition of A. It is the algorithm that gets used behind the scenes in Matlab when you type x=A\b (and A is sparse and square). UMFPACK is free software (GPL)

Also note if y=Ax, and x is known, but y is not, you compute y by performing a sparse matrix vector multiply, not by solving a linear system.

codehippo
  • 1,379
  • 1
  • 8
  • 19
3

Reading the boost documentation, it does not seem like solving w.r.t x is implemented. Solving in y is only a matter of matrix-vector product, which seems implemented in ublas.

One thing to keep in mind is that blas only implement 'easy' operations like addition, multiplication, etc... of vector and matrix types. Anything more advanced (linear problem solving, like your "solve in x y = A x", eigen vectors and co) is part of LAPACK, which built on top of BLAS. I don't know what boost provides in that respect.

David Cournapeau
  • 78,318
  • 8
  • 63
  • 70
3

Boost's linear algebra package's tuning focused on "dense matrices". As far as I know, Boost's package do not have any linear-system-solver. How about use source code in "Numerical Recipe in C (http://www.nr.com/oldverswitcher.html)" ?

Note. There can be subtle index bug in the source code (some code uses array index start from 1)

Jinuk Kim
  • 765
  • 5
  • 5
3

Take a look at JAMA/TNT. I've only used it for non-sparse matrices (you probably want the QR or LU factorizations, both of which have solver utility methods), but it apparently has some facilities for sparse matrices.

Jason S
  • 184,598
  • 164
  • 608
  • 970