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

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


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


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

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

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

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

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

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

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

For generic data the optimal value (in the limit) is zero, which is achieved by any vector with arbitrarily large elements. Just pick a random vector $x_0$ and use $x=\alpha x_0$ where $\alpha$ tends to $\infty$. Am I missing something? The eigenvalue discussion you have is not correct (there is no reason that you should be allowed to add the constraint $y^...


4

Regarding your question about the fastest algorithm for the complete symmetric eigenproblem... There are two main phases to the computation. First is reduction to tridiagonal form using (blocked) Householder reflections. This runs in cubic time w.r.t the matrix dimension, and you can easily code up your own non-blocked version, while the blocked version is ...


4

You'd like to sort eigenvalues/eigenvectors in a way that is continuous as you move through momentum. This highly constrains the sorting - for most k-points, you must sort the eigenvalues to be a close as possible to the eigenvalues from a neighboring k-point. Only at eigenvalue crossings are you allowed to permute the order of bands. Furthermore, in many ...


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


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

Consider what happens when you approximate a function's derivative using finite differences near a boundary. If the boundary is at the point $x_0$, the point $x_1$ is just outside the boundary, then the approximation at the point just before the boundary, $x_{-1}$ is: $$ h^{-4}(f_{-3} -4f_{-2}+6f_{-1}-4f_0+f_1), $$ where $f_0$ and $f_1$ must be given by ...


4

If your matrices are large, why not use a library like ARPACK? The shift-and-invert mode of ARPACK will help you calculate the eigenvalues close to $\sigma$. There are interfaces to ARPACK for most high-level programming languages used for numerical computations (Fortran, C, Python, MATLAB, etc.). For example, a quick tutorial for Python can be found here ...


4

You have noticed that any eigenvector of $A$ is also an eigenvector of your composed matrix. When you compute eigenvectors of your matrix, you can use them in a deflation-type preconditioner, as described in [1, 2]. [1] Frank, J. & Vuik, C. On the construction of deflation-based preconditioners SIAM Journal on Scientific Computing, 2001, 23, 442-462 [2]...


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