15

If your graph is undirected (as I suspect), the matrix is symmetric, and you cannot do anything better than the Lanczsos algorithm (with selective reorthogonalization if necessary for stability). As the full spectrum consists of 100000 numbers, I giess you are mainly interested in the spectral density. To get an approximate spectral density, take the ...


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


13

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


12

Depending on the precision you require for the largest eigenvalue, you could try using the Power Iteration. For your specific example, I would go as far as to not form $A=XX^\mathsf{T}$ explicitly, but compute $x \leftarrow X(X^\mathsf{T}x)$ in each iteration. Computing $A$ would require $\mathcal O(n^3)$ operations whereas the matrix-vector product ...


11

In exact arithmetic you shouldn't need to reorthogonalize regularly, but practically you do. Your u1 and u2 are close to (but not exactly) the true eigenvectors, so your initial deflation almost (but not entirely) removed the true eigenvectors from u3. The tiny components you left behind will be amplified by repeated multiplication by A, you will need to ...


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

In theory, yes. In practice, rounding errors will usually result in (initially slow) convergence to $u_1$. At essentially the same cost one can run the Lanczos algorithm, which will have much faster convergence, and produce the three dominant eigenvalues unless two of these eigenvalues are essentially the same. For Lanczos, selective reorthogonalization is ...


8

According to MathWorld a matrix $A \in \mathbb{R}^{n \times n}$ is positive definite iff $$ (x^T A x) > 0 $$ for all non zero vectors $x\in\mathbb{R}^n$. It is trivial to obtain that $$ x^T\,A\,x = x^T\, \left [ \frac12(A+A^T) \right ]\, x $$ and to recognize that positive definiteness is linked to the spectrum of the symmetric part of matrix $A$. A ...


8

The Lanczos algorithm can be used to put the matrix into tridiagonal form, but it doesn't actually find the eigenvalues and eigenvectors of that tridiagonal matrix. Once you have the matrix in tridiagonal form, the QR algorithm is typically used to find the eigenvalues of the tridiagonal matrix.


7

Calculate the SVD in place of the spectral decomposition. The results are the same in exact arithmetic, as your matrix is symmetric positve definite, but in finite precision arithmetic, you'll get the small eigenvalues with much more accuracy. Edit: See Demmel & Kahan, Accurate Singular Values of Bidiagonal Matrices, SIAM J. Sci. Stat. Comput. 11 (1990)...


7

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

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

Is the matrix also positive, such that eigenvalues are equal to singular values? You can estimate $k$ singular values of a general matrix (possibly nonsquare) with several digits of accuracy provided they are at least somewhat isolated in $8k$ matrix multiplies. The algorithm, as described in the excellent survey by Halko, Martinsson, and Tropp (2012) is: ...


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

A slight generalization of a Householder reflector, as seen in LAPACK's zlarfg, can be used to define a sequence of unitary transformations. In particular, the $j$'th transformation zeros the portion of the matrix below the $j$'th diagonal entry, and also ensures that the $j$'th diagonal entry becomes real. The unblocked algorithm zgeqr2 and the blocked ...


6

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


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

Factor $PB=LU$ yourself and write a routine for evaluating $L^{-1}PAU^{-1}x$ given $x$ (using two backsolves). Then you can solve the problem $(L^{-1}PAU^{-1}-\lambda I)z=0$ with a standard iterative solver for the ordinary eigenvalue problem. If $B$ is singular, compute a left null space basis consisting of the rows of $M$, and a right null space basis ...


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

If only 5 eigenvalues are very significant, the Lanczsos algorithm with $X(X^Tx)$ as matrix-vector multiply should give fast linear convergence after 5 initial steps, hence a fairly accurate largest eigenvalue with few iterations.


5

1) In case of a multiple dominant eigenvalue (and no other of the same absolute value), the power itieration converges to the vector obtained by projecting the starting vector to the dominant eigenspace (if this vector is nonzero). These projections are orthogonal if the matrix is symmetric. Of course, if you start with different starting vectors you'll ...


5

Let $X=\tilde X Z$, where $\tilde X^T D \tilde X=I$. Then we may rewrite $$ (X^T L X)v = \lambda (X^T D X) v $$ as $$ (Z^T \tilde X^T L \tilde X Z) v = \lambda (Z^T Z)v. $$ Now, put $u=Z v$, and premultiply both sides of the previous equation by $Z^{-T}$ to find $$ \tilde X^T L \tilde X u = \lambda u. $$ It is worth noting that, contrary to your ...


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