## Linear equation sys with not positive definit left-hand side

mlimberger
Registered Member
Posts
9
Karma
0
OS

### Linear equation sys with not positive definit left-hand side

Sat Jul 13, 2013 10:42 am
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
ggael
Moderator
Posts
2194
Karma
15
OS

### Re: Linear equation sys with not positive definit left-hand

Sat Jul 13, 2013 11:29 am
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?
mlimberger
Registered Member
Posts
9
Karma
0
OS

### Re: Linear equation sys with not positive definit left-hand

Sat Jul 13, 2013 11:58 am
Hi ggael.

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
ggael
Moderator
Posts
2194
Karma
15
OS

### Re: Linear equation sys with not positive definit left-hand

Sat Jul 13, 2013 8:03 pm
please, try with the devel branch: http://bitbucket.org/eigen/eigen/get/default.tar.bz2. This issue is likely already fixed.
mlimberger
Registered Member
Posts
9
Karma
0
OS

### Re: Linear equation sys with not positive definit left-hand

Sun Jul 14, 2013 11:30 am
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.
ggael
Moderator
Posts
2194
Karma
15
OS

### Re: Linear equation sys with not positive definit left-hand

Sun Jul 14, 2013 1:26 pm
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;}`
mlimberger
Registered Member
Posts
9
Karma
0
OS

### Re: Linear equation sys with not positive definit left-hand

Mon Jul 15, 2013 6:39 am
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.
ggael
Moderator
Posts
2194
Karma
15
OS

### Re: Linear equation sys with not positive definit left-hand

Mon Jul 15, 2013 6:57 am
hm, right SparseLU is expecting a compressed matrix.
mlimberger
Registered Member
Posts
9
Karma
0
OS

### Re: Linear equation sys with not positive definit left-hand

Mon Jul 15, 2013 11:32 am
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.
Dee33
Registered Member
Posts
54
Karma
0
OS

### Re: Linear equation sys with not positive definit left-hand

Mon Jul 15, 2013 2:09 pm
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/
mlimberger
Registered Member
Posts
9
Karma
0
OS

### Re: Linear equation sys with not positive definit left-hand

Mon Jul 15, 2013 3:39 pm
Thank you so much for your support! I opened a Bug session. So let's see...
Dee33
Registered Member
Posts
54
Karma
0
OS

### Re: Linear equation sys with not positive definit left-hand

Tue Jul 16, 2013 2:30 pm
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
tonyhenz
Registered Member
Posts
7
Karma
0

### Re: Linear equation sys with not positive definit left-hand

Fri Jul 19, 2013 8:35 am
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.
ggael
Moderator
Posts
2194
Karma
15
OS

### Re: Linear equation sys with not positive definit left-hand

Fri Jul 19, 2013 3:19 pm
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.
tonyhenz
Registered Member
Posts
7
Karma
0

### Re: Linear equation sys with not positive definit left-hand

Fri Jul 19, 2013 3:39 pm
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.

## Who is online

Registered users: alake, Baidu [Spider], Bing [Bot], boudewijn, colomar, edmael, Exabot [Bot], garthecho, Google [Bot], google01103, Hans, jitseniesen, ken300, koriun, mayankasthana, mgraesslin, mixoff, mmistretta, pedrorodriguez, scummos, SeaJey, SecretCode, Sentynel, Sogou [Bot], urgo, Yahoo [Bot]