# Tag Info

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 ...

6

No, in the general case, there is no suitable workaround. C++ is a statically typed language, and the compiler needs to know all types at compilation time. If your code worked, the following would also const size_t n = std::rand(); Eigen::Matrix<double, n, n> A; and rand() only gives a random number at run-time, which, at compile-time is unknown. ...

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

5

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

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 ...

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

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

Certainly. There are a few things you have to define for your type that are listed on this page in the documentation: https://eigen.tuxfamily.org/dox/TopicCustomizing_CustomScalar.html It basically boils down to defining arithmetic operators appropriately for your type, plus specializing a traits template NumTraits that describes your type. The link above ...

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 ...

3

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 ...

3

This problem is known as joint diagonalization, and it has two variants: orthogonal, in which the basis vectors are orthonormal, and non-orthogonal, which is harder to solve, but which may be more appropriate to your application. The simplest method I know of seeks a unitary matrix $U$ that minimizes the sum of squares of the off-diagonal elements of $U^HAU$...

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

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 ...

2

As its documentation suggests in multiple places, rref is "mainly of academic interest" (read: "used only to explain Gaussian elimination to undergrads"), and is not a serious competitor of SVD-based algorithms in terms of stability. I recommend against it.

1

I get it. template <typename Number> // e.g. double or complex Eigen::Matrix<Number, Eigen::Dynamic, Eigen::Dynamic> image_COD( const Eigen::Matrix<Number, Eigen::Dynamic, Eigen::Dynamic>& M) { Eigen::CompleteOrthogonalDecomposition< Eigen::Matrix<Number, Eigen::Dynamic, Eigen::Dynamic>> cod(M); const Eigen::...

1

You can use arpack , 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  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 ...

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