14

The convergence behavior you are seeing is actually expected. One of things that makes the Lanczos method so interesting is that it does a good job of simultaneously converging eigenvalues at both ends of the spectrum. I assume your expectation of converging only the largest eigenvalues is based on the fact that, as expected from the Power iteration ...


13

Compute the largest-magnitude eigenvalue $\lambda_{max}$ of $A$ (with, say, eigs('lm')). Then compute the largest magnitude (negative) eigenvalue $\hat{\lambda}_{max}$ of $M = A - \lambda_{max}I$ (again, through a standard call to eigs('lm')). Observe that $\hat{\lambda}_{max} + \lambda_{\max} = \lambda_{min}(A)$. The reason why this holds is explained here. ...


11

For eigenvalues, simply take $k$ largest or smallest eigenvalues of $T$. They are good approximations of $A$, provided that the number of Lanczos iterations is large compared to $k$. Things are a little trickier if we want eigenvectors as well. The simplest way is to multiply each eigenvector $\mathbf{u}_i$ of $T$ by $V$ to the left, where $V$ is, as you ...


10

It's not implementation-dependent in the sense that this is a mathematical operation performed on your matrix. However, it is very much matrix-dependent. If your matrix is diagonalizable and $A=XDX^{-1}$, then zeroing out some element adds a small perturbation matrix $E$, so the new eigenvalues will be (assuming the matrix $X$ does not change much) $$ X^{-1}...


9

You should specify the eigenvalues you want with which="SM", for example. Check the following snippet. I also changed the solver, since your system is symmetric. import numpy as np from scipy.sparse.linalg import eigsh import matplotlib.pyplot as plt n = 200 h = 2/(n-1) # domain for x and y is [-1, 1] L = np.diag(np.ones(n-1), k=-1) - np.diag(2*np....


8

Your main concern is not destroying sparsity - a transformation does not necessarily destroy it. Consider the transformed eigenvalue problem, that eliminates the constraint $x^Tx = 1$ $$ \min_x \frac{x^TAx}{x^Tx} \qquad s.t.\quad Bx = 0$$ Consider a matrix $Z$ that is a basis for the nullspace of $B$, i.e. $BZ = 0$. Then, restricting $x$ to be in the ...


7

FILTLAN is a C++ library for computing interior eigenvalues of sparse symmetric matrices. The fact that there is a whole package devoted to just this should tell you that it's a pretty hard problem. Finding the largest or smallest few eigenvalues of a symmetric matrix can be done by shifting/inverting and using the Lanczos algorithm, but the middle of the ...


7

The following paper suggests that the Jacobi-Davidson method can be used to target eigenvectors based on "any property that can be computed from the eigenvector", which would seem to include overlap with a given vector. https://journals.aps.org/prb/abstract/10.1103/PhysRevB.66.245104 According to the paper, the key is just to reorder the QR decomposition ...


6

Arnold Neumaier's answer is discussed in more detail in section 3.2 of the paper "Approximating Spectral Densities of Large Matrices" by Lin Lin, Yousef Saad and Chao Yang (2016). Some other methods are also discussed but the numerical analysis at the end of the paper shows that the Lanczos method outperforms these alternatives.


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

A 100,0000 by 100,000 symmetric dense matrix in single precision requires 20 gigabytes of memory (storing only the upper triangle) or 40 gigabytes of memory for double precision. Thus it is too large to fit within the memory of available GPU's. In order to solve this problem using GPU acceleration you'd have to develop an algorithm that sends smaller ...


6

Yes, for matrices A, B, the last column of matrix product A*B can be written as A*(the last column of B). You can use this fact to get the last column of the reconstruction using only O(N^2).


6

There is a field of study known as eigenvalue sensitivity analysis or eigenvalue perturbation analysis that allows you to estimate the effect of small matrix perturbations on the eigenvalues and eigenvectors. The basic technique used for this is differentiating the eigenvalue matrix equation, $$AX = X\Lambda.$$ For situations where the eigenvalues of the ...


6

This is something like a truncated SVD or eigenvector expansion of your solution. If you take $$x_m = \sum_{j=1}^m \frac{q_j\cdot b}{\lambda_j}q_j$$ with $m=n$, this is the exact solution to $Ax=b$. If you take instead the $m<n$ eigenvalues in the sum instead, you get an approximation. The problem is that it's a pretty poor approximation; since your ...


6

I don't know for sure whether there is an existing name for this method, but @jessechan's suggestion of "truncated eigenvector expansion" sounds perfectly fine to me (and most people would understand it). Your question had a second part embedded, namely why is this not what everyone does to get a good initial guess for iterative solvers? The answer to this ...


6

For a symmetric 3x3 matrix, one Householder transformation will bring your matrix in tridiagonal form. The required algorithm is given (for general $n\times n$ matrices) on page 459 of Matrix Computations, 4th edition, Algorithm 8.3.1. For a $3\times 3$ matrix, it's just one Householder reduction instead of a loop. For the subsequent tridiagonal matrix, ...


5

10000 by 10000 isn't really very big. A double precision matrix of this size requires 800 megabytes of memory. You'll need at least as much memory to hold the matrix of resulting eigenvectors, and you'll need additional working storage, but this is well within the capability of a typical desktop machine with 8 gigabytes or more of memory. It would also ...


5

With the reverse communication protocol of ARPACK, you do not need to store the $3n \times 3n$ matrix explicitly: you just need to provide two functions that compute: $ \left[ \begin{array}{c} x \\ y \\ z \end{array}\right] \rightarrow \left[ \begin{array}{c} -A_0 x \\ y \\ z \end{array}\right]$ and $ \left[ \begin{array}{c} x \\ y \\ z \end{array}\right] \...


5

If, as you say, you are sure that you have a symmetric-definite pencil (that is, $\mathbf A$ is symmetric, and $\mathbf B$ is symmetric positive-definite), then LAPACK already has something for directly handling your problem: dsygv(). What it does is to perform a Cholesky decomposition of $\mathbf B$ (if in fact your $\mathbf B$ is not symmetric positive-...


5

There has been some good research on this recently. The new approaches use "randomized algorithms" which only require a few reads of your matrix to get good accuracy on the largest eigenvalues. This is in contrast to power iterations which require several matrix-vector multiplications to reach high accuracy. You can read more about the new research here: ...


5

Yes. Simply run $n$ steps of column-pivoted QR on the adjoint of your matrix. You may want to look into zgeqp3. If you express the result of the pivoted $QR$ decomposition as: $$ A^H \Omega = Q R, $$ your $B$ matrix will be the adjoint of the left-most $n$ columns of $A^H \Omega$, as $\Omega$ is a permutation matrix.


5

Let $w_k$ be the k-th column of W (the kth eigen-vector) and $v_k$ be the k-th element of v (the kth eigen-value). Then we can write: $$ G = \sum_k w_k v_k w_k^T $$ which, element-wise, is equivalent to: $$ G_{i,j} = \sum_k v_k w_{i,k} w_{j,k} $$ where the indices are column-major. Since we always get $w_k$ twice, an overall sign on a column of $W$ cancels. ...


5

MATLAB always uses the LAPACK libraries to calculate eigenvectors which works on double precision, floating point numbers - MATLAB's default data type. Mathematica's method depends on its input type. For example, when you do TestMatrix = {{1, 2, 3}, {3, 1, 2}, {2, 3, 1}} Eigenvectors[TestMatrix] You'll get an exact answer involving Sqrt[3] and so on. ...


5

The algorithm used by eigs is called Arnoldi iteration (in shift-and-invert mode). You could think of it as a sophisticated way of doing an inverse power iteration. The inverse power iteration method works as follows: if $A$ is your matrix and you know that the Perron root (the spectral radius of $A$) is $\lambda$, then a power iteration on the matrix $B = ...


5

I suspect the main problem is the magnitude of the values. If you divide through by $\cosh$ to make all the numbers smaller, then ApproxFun doesn't seem to have a problem finding all the roots. function lhs(ε, α) A, B = sqrt((α-1)/(2ε)), sqrt((α+1)/(2ε)) return sqrt(α*α-1)*(cos(A) - 1/cosh(B)) - sin(A)*tanh(B) end u = Fun(x -> lhs(1e-3, x), [1.0,...


4

From geometrical tools (www.geometrictools.com) http://www.geometrictools.com/Documentation/EigenSymmetric3x3.pdf I hope this helps


4

A Google search returns many relevant links for this problem since it is a common problems in mechanics. A simple algorithm for symmetric and real 3x3 matrices is presented on Wikipedia: http://en.wikipedia.org/wiki/Eigenvalue_algorithm#3.C3.973_matrices


4

MATLAB's eig() function normalizes the eigenvectors to have unit norm; it doesn't normalize the eigenvalues. You can check the norm of the vectors to verify their correctness. Remember, if $\bf x$ is an eigenvector of $A$, then so is $c\,\bf x$. You can "turn off" the normalization by using [v,d] = eig(A,'nobalance') but this also disables ...


4

If you are only interested in the smallest eigenvalue, the conjugate gradient method applied to the matrix $L$ gives you a good approximation after a reasonably small number of steps, and you won't have to solve any linear systems. The details are in Y. Saad's book on iterative methods, but here a short summary: From the coefficients that are computed in ...


4

Assuming your overall matrix is positive definite (it is definitely symmetric), then I would suggest looking into algebraic multigrid (AMG) methods as preconditioners. They compute hierarchies of sparsified matrices themselves. If you're already using PETSc, take a look at the hypre preconditioner. Using this may force you to actually multiply out the ...


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