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


14

It depends a lot on the size of your matrix, in the large-scale case also on whether it is sparse, and on the accuracy you want to achieve. If your matrix is too large to allow a single factorization, and you need high accuracy, the Lanczsos algorithm is probably the fastest way. In the nonsymmetric case, the Arnoldi algorithm is needed, which is ...


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


12

In gerschgorin's theorem, the diagonal entries $A_{ii}$ of the matrix are the eigenvalue estimates, and the radii $r_i$ of the Gerschgorin disks are corresponding error bounds. Thus $\min_i A_{ii}-r_i$ is a lower bound on the eigenvalues, and $\max_i A_{ii}+r_i$ is an upper bound. Note that these bounds are generally poor unless the off-diagonal entries are ...


12

The fastest method will likely depend upon the spectrum and normality of your matrix, but in all cases Krylov algorithms should be strictly better than power iteration. G.W. Stewart has a nice discussion of this issue in Chapter 4, Section 3 of Matrix Algorithms, Volume II: Eigensystems: The power method is based on the observation that if $A$ has a ...


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

We have the matrix Laplacian matrix $G=A^TA$ which has a set of eigenvalues $\lambda_0\leq\lambda_1\leq\ldots\leq \lambda_n$ for $G\in\mathbb{R}^{n\times n}$ where we always know $\lambda_0 = 0$. Thus the Laplacian matrix is always symmetric positive semi-definite. Because the matrix $G$ is not symmetric positive definite we have to be careful when we ...


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


10

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


9

TQL cannot be parallelized. The standard parallel algorithm is that of Cuppen: JJM Cuppen, A divide and conquer method for the symmetric tridiagonal eigenproblem, 1980. http://www.springerlink.com/content/t21365q2gh702714/ see also: F. Tisseur, A parallel divide and conquer algorithm for the symmetric eigenvalue problem on distributed memory ...


9

I think that you are confusing the unique matrix square-root of Hermitian positive semi-definite matrix $A$, i.e., a Hermitian positive semi-definite matrix $B$ satisfying, $$ B^2 = A, $$ with the non-unique problem of finding a matrix $C$ satisfying $$ C^H C = A, $$ where clearly the mapping $C \mapsto Q C$, for any unitary $Q$, preserves the ...


8

Youssef Saad's book Numerical Methods for Large Eigenvalue Problems, 2nd edition uses the norm of the residual vector to define convergence criteria. He defines the residual vector as follows on page 59: Given a matrix $\mathbf{A} \in \mathbb{C}^{n \times n}$, a putative eigenvalue $\widetilde{\lambda} \in \mathbb{C}$ and a putative eigenvector $\widetilde{\...


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

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

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.


8

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 Power Iteration (or Power Method), e.g. what Dan is describing, should always converge, albeit at the rate $\left|\lambda_{n-1}/\lambda_{n}\right|$. If $\lambda_{n-1}$ is close to $\lambda_n$, it will be slow, but you can use extrapolation to get around that. It may seem complicated, but an implementation in pseudo-code is given in the paper.


7

Gershgorin's Theorem provides a bound on where to find each eigenvalue, but it doesn't provide an algorithm to actually calculate them. You can take an initial guess $x_0$ within each Gershgorin disk, and use the shifted inverse power method to find the eigenvalue closest to $x_{0}$. But there is no guarantee that using one initial guess per Gershgorin ...


7

For simplicity, assume that there is only one parameter $t$ rather than your two. In order that you can have continuous eigenspaces, you need to assume that the associated eigenvalues do not nearly cross. (For nearly crossing eigenvalues it may very well happen that the eigenspaces are essentially exchanged though the eigenvalue curves do not touch. This ...


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

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

In essence, what you're asking is to find the square root A of a matrix G, so that $$G = A^T A.$$ There are many ways to do that if $G$ is a symmetric matrix. For example, if $G$ is symmetric, then the Cholesky decomposition $G=L^TL$ provides you with one answer: $A=L$. But you already found another answer, with the matrix $A$ you already have. What this ...


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

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

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


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