Hot answers tagged

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

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);


5

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


5

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


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

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


4

Despite my initial conviction that I was including all relevant compiler flags, in fact I wasn't. I forgot the -fopenmp flag on the C++ compiler, which enables multithreading and a fairer comparison between implementations. I also found that -march=native makes the code three times faster than the more general -march=nocona which I was originally using, and ...


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

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


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


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

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

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

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


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