Hot answers tagged

12

The LU factors of a sparse matrix are at least somewhat sparse. The $Q$ matrix in QR can also somewhat preserve sparsity, and is typically used when the matrix is very long and skinny. The SVD of a sparse matrix will almost always have fully dense $U$ and $V$ factors, so it destroys any reason to perform the computations treating the matrix sparsely.


8

There are two relatively convenient options for calculating selected (e.g. a few largest or smallest) eigenvalues using Eigen. The first is Spectra, a header-only C++ library based on Eigen that uses algorithms similar to ARPACK (implicitly-restarted Arnoldi) to calculate a few eigensolutions. Since it is header-only, you simply download and include the ...


7

You have a size mismatch issue: A is a count x 2*count matrix, and you are trying to solve Ax=B with B a 2*count x 1 vector. Moreover, if you compile without -DNDEBUG, you should get a nice assertion telling you this is wrong. To resize a matrix or vector: A.resize(count, 2*count); B.resize(count);


6

You need to declare a LDLT object like this: LDLT<MatrixXd> ldlt(A); x = ldlt.solve(b); y = ldlt.solve(x); ...


6

You're looking for the SelfAdjointEigenSolver class, and there is also an example in the user manual that I report here: #include <iostream> #include <Eigen/Dense> using namespace std; using namespace Eigen; int main() { Matrix2f A; A << 1, 2, 2, 3; cout << "Here is the matrix A:\n" << A << endl; ...


6

We have the matrix $A$ that can be expressed as $A = JGJ^T$. The first thing is to calculate the QR decomposition of matrix $J$. Because of the low rank of the matrix it can be done very fast with, for instance, modified Gram Schmidt algorithm. Now we can write $A$ as $A = QR G R^TQ^T$, where $Q$ is an orthonormal matrix ($Q^T Q = I$). We define $F$ as ...


5

The other answers already tell you what went wrong, but I will add a terminology note: the term for what is happening is that the pencil $A - \lambda B$ is a singular matrix pencil, i.e., $\det (A - \lambda B)$ is identically equal to zero. So there are no generalized eigenvalues (or, at least, they cannot be defined as usual as the roots of the generalized ...


5

Please make following changes in your code and check the answer b(i1,0)=1; b(i1,1)=a(i1); b(i1,2)=a(i1)*a(i1)+1.0/12.0; b(i1,3)=pow(a(i1),3)+a(i1)/4; b(i1,4)=1.0/80.0+pow(a(i1),2)/2+pow(a(i1),4);


4

This is a O(m) problem where m is the size of the big array. So this should be fast: int n = v[v.size()-1]; int k=0; VectorXi lengths[n+1], starts[n+1]; for(int i=0; i<v.size(); ++i) { int i0 = i; while(v[i]==k && i<n) ++i; lengths[k] = i-i0; starts[k] = i; ++k; }


4

As documented here: x = A.triangularView<Upper>().solve(b); or x = A.triangularView<Lower>().solve(b);


4

This has nothing to do with Eigen: if you cannot, just by looking at a matrix, determine how eigenvalues should be labelled, you cannot expect Eigen to do it for you. Also, it is up to you to define precisely, mathematically, how the labels should be assigned. Judging from your example plot, you seem to be assuming that labelled eigenvalues form smooth ...


4

I've rolled my own. Here is a MCVE: #include <Eigen/Core> #include <Eigen/Sparse> #include <iostream> #include <fstream> #include <vector> using namespace Eigen; typedef Triplet<int> Trip; template <typename T, int whatever, typename IND> void Serialize(SparseMatrix<T, whatever, IND>& m) { std::...


4

Eigen's LDLT class actually performs a U^T.D.U factorization if the input matrix is column-major with the symmetric entries stored in the upper triangular part: MatrixXd A; // fill at least the upper triangular part of A LDLT<MatrixXd,Upper> udu(A); // Only the upper part of A is read to form U udu.matrixU(); // returns an expression of the triangular ...


3

As others have mentioned in comments, solving $A^TAx=A^Tb$ (the so-called least-norm solution) can be performed without explicitly forming $A^TA$ with the QR decomposition. It works like this: $$A=QR$$ where $Q$ has the interesting properties that $Q^TQ=I$ and $R$ is square and triangular. Using this information, the least squares system $A^TAx=A^Tb$ is ...


3

I suspect the root of your trouble is what has been detected in the comments by Vibe: For any number $\omega\in \mathbb{K}$ (with $\mathbb{K}= \mathbb{R}$ or $\mathbb{C}$) you can find $\boldsymbol{X}$ such that $AX = \omega BX$ (with $A$ and $B$ taken in your concrete example). You have already decomposed the problem in 4 blocks of 3 variables. Then let us ...


2

This is a typical case where it does not make sense to actually store the matrix. Rather, consider the matrix to be an operator where if you need to multiply by it, i.e., form $y=Ax$, you do $y=b (b\cdot x)$ instead if you need to access an element $A_{ij}$, you instead compute it as $b_ib_j$. In other words, knowing what $A$ is, you should just not store ...


2

Given that the accepted answer has typos and does not give full code, I decided to post mine. #include <iostream> #include <Eigen/Dense> using namespace std; using namespace Eigen; int main() { MatrixXd A; A.resize(3,3); VectorXd b; b.resize(3); A << 13, 5, 7 , 5 , 9, 3 , 7 , 3, 11; b << 3, ...


2

Eigen::SparseLU<Eigen::SparseMatrix<double> > solverA; A.makeCompressed(); solverA.analyzePattern(A); solverA.factorize(A); if(solverA.info()!=Eigen::Success) { std::cout << "Oh: Very bad" <<"\n"; } else{ std::cout<<"okay computed"<<"\n"; } Eigen::VectorXd solnew = solverA.solve(b); You don't need to mention ...


2

Based on your comments, the most likely explanation is as follows. In dggevx.f, there is the following paragraph: * Optionally also, it computes a balancing transformation to improve * the conditioning of the eigenvalues and eigenvectors (ILO, IHI, * LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for * the eigenvalues (RCONDE), and ...


1

You can use arpack [1], that implements the Arnoldi algorithm for computing eigenpairs. Since arpack communicates with client code only through matrix vector product, you can use your own matrix type (including eigen sparse matrix). There is also a version with C++ bindings [2] that may be easier to use with eigen. To compute the smallest eigenvalue, it may ...


1

If one of the matrices is invertible, you can convert it to a standard eigenvalue problem $(B^{-1}A) x = \lambda x$. With Eigen, you can thus do: EigenSolver<MatrixXcd> eig(B.lu().solve(A)); and, if needed, you can also re-normalize the eigenvectors so that $x^T B x = 1$.


1

EDIT: I just noticed a glitch in your code. You used selfadjoint and triangular. Try this: L.triangularView<Lower>().transpose().solve(b); The rest of my answer may still be useful, so I won't delete it. Presumably, $L$ uses the column-major storage order. If not, just swap row for column everywhere in this post. $L^T$ can be represented by casting $L$...


1

I would second the opinion expressed in the context. For that small problem and limited usage, you don't need a library. Generation of a structured grid and retrieving points can be coded up in at most several screens of code. However, I would point out for you ViennaGrid library, which can be used and actually provides STL iterators. In addition, the ...


1

The answer to your question is benchmarking. Both Armadillo and Eigen provide some benchmarks in their documentation. The only problem is that they don't compare to the same libraries, but you can still get an idea. Maybe there are more comprehensive benchmarks elsewhere... In any case, the overhead of wrapping LAPACK in a C++ library will be much smaller ...


1

Numerical computation of Generalized Complex Schur decomposition can be performed by calling zgges() LAPACK function. For example, see NETLIB zgees documentation, or a documentation for any other BLAS/LAPACK library implementation. Eigen is technically nothing else, but a very convenient templated library of wrappers and algorithms, also including some ...


1

After your second edit, you are already here: for (int i = 0; i < Ax.size(); i++) { double Ax_i = 0.0; for (int dataIdx = Aindptr[i]; dataIdx < Aindptr[i + 1]; dataIdx++) { Ax_i += Adata[dataIdx] * x[Aindices[dataIdx]]; } Ax[i] = Ax_i; } But notice how you in row i+1 you are initializing dataIdx to the value Aindptr[i+1] it already ...


1

I do not have too much experience using eigen, but when you are solving the systems with the SVD decomposition, actually you are doing the followin: $ A x = USV^T x = b $ and you use the SVD decomposition of A to isolate x $x = VS^{-1}U^T b$. The same for the other system: $ A^T y = VSU^T y = c $ and you use the SVD decomposition of A to isolate x $y =...


1

An LDLT decomposition will give you the information you need. If the diagonal $D$ matrix from the LDLT decomposition of your matrix is positive definite, your matrix is also positive definite. Here is a link to some documentation discussing the LDLT decomposition in Eigen: http://eigen.tuxfamily.org/dox-2.0/TutorialAdvancedLinearAlgebra.html


1

This is an old question and it does not contain enough information to answer it with certainty. But since I already stumbled on it once, I want to propose a list of actions which should be able to help resolve it with all the uncertainties. The very first action to do is to check if the original matrix itself contains valid entries. hasNaN() and allFinite() ...


Only top voted, non community-wiki answers of a minimum length are eligible