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


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

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

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

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

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

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

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


5

Yes. If $A\in\mathbb{R}^{n\times n}$ is symmetric, there exists an orthogonal basis of eigenvectors $v_1,\dots,v_n$ (sorted by decreasing magnitude) of $A$. If you expand your starting vector $u_3$ with respect to this basis, you get $$ u_3 = \alpha_1 v_1 + \alpha_2 v_2 + \dots + \alpha_n v_n,$$ where $\alpha_1 = \alpha_2 =0$ because $u_3$ is orthogonal to $...


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

When $B$ is indefinite, the eigenvalues may be complex, and there is little advantage exploiting the symmetry. Instead, one generally uses the QZ algorithm. Edit: If you know a priori that all your eigenvalues are real, it is probably because you can establish a priori that some linear combination $C=sA+tB$ is positive definite. In this case, you should ...


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

Eigenvector solutions to the standard eigenvalue problem formulation ($A\mathbf{v}=\lambda\mathbf{v}$) are unique up to a constant multiplier. It is clear that if there exists an eigenvector solution $\mathbf{v_a}$ which satisfies $A\mathbf{v_a}=\lambda\mathbf{v_a}$, then for any real\complex scalar $k$, $\mathbf{v_b}=k\mathbf{v_a}$ is also valid simply ...


5

The computation of integer-valued matrix determinants has been a subject of considerable research. Using exact arithmetic the Smith normal form can be computed, and from this diagonal form the determinant is easily found. Saunders and Wan (2004), Smith Normal Form of Dense Integer Matrices, Fast Algorithms into Practice, say "Over the past thirty years, ...


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

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


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