Hello.

So far, I used SimplicialLDLT to solve for linear equations system with a sparse, symmetric, positive definite left-hand side (LHS).

Now, the LHS has to be extended and becomes sparse, symmetric but not positive definite anymore. The size is around 15000 x15000.
The structure of LHS has the characteristic, that all diagonal elements but one have non-zero entries: diag(LHS) = [x,x,x,x,x,x,x,x,...,x,x,0]
The RHS then is [x,x,x,x,x,...,x,x,0].

Therefore I switched from SimplicialLDLT (which is only valid for symmetric positive definite problems, right?) to SparseLu from
which I thought it is mybe the best way for solving. Can you agree or do you think there are better alternatives?

Based on x = LHS^(-1)*RHS, my implementation now looks like:

-----------
...
#include "Eigen/SparseLU"
#include "Eigen/OrderingMethods"
...

x = Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int> >(LHS).solve(RHS);

----------

Unfortunetaly I receive the following error message:

Eigen::DenseCoeffsBase<Derived, 1>::Scalar& Eigen::DenseCoeffsBase<Derived, 1>::operator[](Eigen::DenseCoeffsBase<Derived, 1>::Index) [with Derived = Eigen::Ref<Eigen::Matrix<int, -0x00000000000000001, 1>, 0, Eigen::InnerStride<1> >, Eigen::DenseCoeffsBase<Derived, 1>::Scalar = int, Eigen::DenseCoeffsBase<Derived, 1>::Index = long int]: Assertion `index >= 0 && index < size()' failed.

Sounds for me like the reson for the error is based on the COLAMDOrdering. I also testet NaturalOrdering and AMDOrdering, but the same error message appeared.

I would greatly appreciate it if you could give me some help, support or comments on how to solve this kind of least squares in an efficient way.

Thank you very much in advance!
Marco

Do you use the devel branch? or the 3.2-beta-1? There has been a recent fix in SparseLU. If that still fails with the latest revision, could you send the matrix?

Hi ggael.

Thanks for your fast reply! I use the Eigen 3.2-beta-1 version.

I tested SparseLu with a small 3x3 LHS test matrix where the same problem appears:

Eigen::MatrixXd lhs(3,3);
lhs.coeffRef(0,0)=-3; lhs.coeffRef(0,1)=11; lhs.coeffRef(0,2)=-2;
lhs.coeffRef(1,0)=11; lhs.coeffRef(1,1)=-5; lhs.coeffRef(1,2)=-2;
lhs.coeffRef(2,0)=-2; lhs.coeffRef(2,1)=-2; lhs.coeffRef(2,2)=0;

and 3x1 RHS:
Eigen::VectorXd rhs(3,1);
rhs<<1,1,1;


In case of using a dense matrix, the result yields:

x = lhs.inverse()*rhs =
[-0.266667, -0.233333, -1.38333]

---

But when I use SparseLu with

Eigen::SparseMatrix<double>lhs(3,3);
... same entries...

and

x = Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int> >(lhs).solve(rhs);

the error message appears.

Any idea why SparseLu cannot solve it?
Or do you anyway think that SparseLu is not the the best way for such a linear equations system?

Regards, Marco

please, try with the devel branch: http://bitbucket.org/eigen/eigen/get/default.tar.bz2. This issue is likely already fixed.

Same Error with the devel branch Hmmm...Do you have any suggestion what other solver I could try?
The important thing is that it works properly also for relatively large systems.

I can reproduce the issue with the beta-1 but it does work fine with the devel branch so make sure you properly configured your compiler to use the devel branch and not another copy of Eigen.

Here is my test program:

Code: Select all

#include <Eigen/Dense>
#include <Eigen/SparseLU>
using namespace Eigen;

int main()
{
  MatrixXd lhs(3,3);
  lhs.coeffRef(0,0)=-3; lhs.coeffRef(0,1)=11; lhs.coeffRef(0,2)=-2;
  lhs.coeffRef(1,0)=11; lhs.coeffRef(1,1)=-5; lhs.coeffRef(1,2)=-2;
  lhs.coeffRef(2,0)=-2; lhs.coeffRef(2,1)=-2; lhs.coeffRef(2,2)=0;

  VectorXd rhs(3), x;
  rhs << 1,1,1;

  x = lhs.inverse()*rhs;
 
  std::cout << "dense:  " << x.transpose() << " -> " << (lhs*x-rhs).norm() << "\n";

  SparseMatrix<double> slhs(lhs.sparseView());
  SparseLU<SparseMatrix<double>, COLAMDOrdering<int> > slu(slhs);
  x = slu.solve(rhs);
 
  std::cout << "sparse: " << x.transpose() << " -> " << (slhs*x-rhs).norm() << "\n";
 
  return 0;
}


The compiler was configured correctly. But you are right, your Testprogram is working fine. My testprogram was:

Code: Select all

Eigen::SparseMatrix<double> lhs(3,3);
lhs.coeffRef(0,0)=-3;lhs.coeffRef(0,1)=11;lhs.coeffRef(0,2)=-2;
lhs.coeffRef(1,0)=11;lhs.coeffRef(1,1)=-5;lhs.coeffRef(1,2)=-2;
lhs.coeffRef(2,0)=-2;lhs.coeffRef(2,1)=-2;lhs.coeffRef(2,2)=0;

Eigen::VectorXd res; rhs<<1,1,1;

Eigen::VectorXd res;
Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int> > slu(lhs2);
res = slu.solve(rhs);


In this case, the error message appears. But I realized that if I add an additional lhs.makeCompressed(); is works:

Code: Select all

Eigen::SparseMatrix<double> lhs(3,3);
lhs.coeffRef(0,0)=-3;lhs.coeffRef(0,1)=11;lhs.coeffRef(0,2)=-2;
lhs.coeffRef(1,0)=11;lhs.coeffRef(1,1)=-5;lhs.coeffRef(1,2)=-2;
lhs.coeffRef(2,0)=-2;lhs.coeffRef(2,1)=-2;lhs.coeffRef(2,2)=0;

lhs.makeCompressed();

Eigen::VectorXd res; rhs<<1,1,1;

Eigen::VectorXd res;
Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int> > slu(lhs2);
res = slu.solve(rhs);


For the real case with the large matrix, I initialize the matrix also with the expected size and add single elements by LHS.coeffRef(i,j).
Can this cause an error when using SparseLu? At least, the error still remains here, even if I compress the matrix.

hm, right SparseLU is expecting a compressed matrix.

Unfortunately, it still does not work for the larger matrix. The error message that arises is:

const Eigen::internal::solve_retval<Eigen::SparseLU<_MatrixType, _OrderingType>, Rhs> Eigen::SparseLU<_MatrixType, _OrderingType>::solve(const Eigen::MatrixBase<OtherDerived>&) const [with Rhs = Eigen::Matrix<double, -0x00000000000000001, 1>, _MatrixType = Eigen::SparseMatrix<double>, _OrderingType = Eigen::COLAMDOrdering<int>]: Assertion `m_factorizationIsOk && "SparseLU is not initialized."' failed.

Using dense matrices with

x = LHS.jacobiSvd(Eigen::ComputeThinU|Eigen::ComputeThinV).solve(RHS);

works.

Is there a convenient way that I send you the LHS matrix (733x733) and RHS vector (733x1)? I am really stucked with this problem.

Please use the functions saveMarket() and saveMarketVector() defined in unsupported/Eigen/src/SparseExtra/MarketIO.h
Then file a bug here and attach the required files : http://eigen.tuxfamily.org/bz/

Thank you so much for your support! I opened a Bug session. So let's see...

ggael wrote:hm, right SparseLU is expecting a compressed matrix.

SparseLU now accepts non compressed matrices :
https://bitbucket.org/eigen/eigen/commits/02d91f29a8ff/
Changeset: 02d91f29a8ff
User: dnuentsa
Date: 2013-07-16 15:15:53
Summary: Fix Sparse LU for matrices in non compressed mode
Affected #: 1 file

For following posts : visit http://eigen.tuxfamily.org/bz/show_bug.cgi?id=630

Hi Dee33 and all,
Now i have met exact the same problem. According to the previous discussions, i have used the newest version of Eigen, for example 3.2rc-1. However, this issue still exists. But it is a little strange, i have tested my problem, if the sparse matrix dimension is small, to say 70000 * 70000, the SparseLU works fine, but when the dimension increases to 200000 * 200000, the error presented above will show again. I have tested my problem by using BiCGSTAB solver, it works fine but quite slowly. Any suggestions will be appreciated.

You mean you get this exact same error message:

Code: Select all

Eigen::DenseCoeffsBase<Derived, 1>::Scalar& Eigen::DenseCoeffsBase<Derived, 1>::operator[](Eigen::DenseCoeffsBase<Derived, 1>::Index) [with Derived = Eigen::Ref<Eigen::Matrix<int, -0x00000000000000001, 1>, 0, Eigen::InnerStride<1> >, Eigen::DenseCoeffsBase<Derived, 1>::Scalar = int, Eigen::DenseCoeffsBase<Derived, 1>::Index = long int]: Assertion `index >= 0 && index < size()' failed.


?

You can also check you are really using 3.2-rc1 by printing cout << EIGEN_MAJOR_VERSION << "." << EIGEN_MINOR_VERSION; You should get 1.92.

ggael wrote:You mean you get this exact same error message:

Code: Select all

Eigen::DenseCoeffsBase<Derived, 1>::Scalar& Eigen::DenseCoeffsBase<Derived, 1>::operator[](Eigen::DenseCoeffsBase<Derived, 1>::Index) [with Derived = Eigen::Ref<Eigen::Matrix<int, -0x00000000000000001, 1>, 0, Eigen::InnerStride<1> >, Eigen::DenseCoeffsBase<Derived, 1>::Scalar = int, Eigen::DenseCoeffsBase<Derived, 1>::Index = long int]: Assertion `index >= 0 && index < size()' failed.


?

You can also check you are really using 3.2-rc1 by printing cout << EIGEN_MAJOR_VERSION << "." << EIGEN_MINOR_VERSION; You should get 1.92.

Hi ggael,
Thanks for your prompt reply. I have printed the output, which is 1.92. And i found the error message is not exactly the same, but slightly different, please see below. The strange thing is, for a different matrix dimension it works quite nice, so i am confused.

../eigen3/Eigen/src/Core/PlainObjectBase.h:265: void Eigen::PlainObjectBase<Derived>::resize(typename Eigen::internal::traits<T>:$$ts<T>::Index) [with Derived = Eigen::Matrix<int, -0x00000000000000001, 1, 0, -0x00000000000000001, 1>]: Assertion `((SizeAtCompileTime == D$$e == Dynamic && (MaxSizeAtCompileTime==Dynamic || size<=MaxSizeAtCompileTime)) || SizeAtCompileTime == size) && size>=0' failed.