20

What's your working definition of "positive semidefinite" or "positive definite"? In floating point arithmetic, you'll have to specify some kind of tolerance for this. You could define this in terms of the computed eigenvalues of the matrix. However, you should first notice that the computed eigenvalues of a matrix scale linearly with the matrix, so that ...


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

In a nutshell, the QR algorithm applied to a matrix $A$ is an iterative procedure that converges to the real Schur decomposition: a unitary matrix $Q$ and a matrix $R$ in block upper triangular form (see below) such that $A = QRQ^T$. It follows that the columns of $Q$ are the eigenvectors (which are the principal objects that are computed!) and that $R$ has ...


12

The degeneracy of some eigenvalues looks to me like the hallmark of the breakdown of the Lanczos algorithm. The Lanczos algorithm is one of the more commonly used methods to approximate the eigenvalues and eigenvectors of Hermitian matrices; it's what scipy.eigsh() uses, through a call to the ARPACK library. In exact arithmetic, the Lanczos algorithm ...


12

Finding the eigenvalues for the Schrödinger equation is really similar to finding the eigenvalues for the wave equation. You start with your differential equation $$\left[-\frac{1}{2}\nabla'^2 + V(r)\right]\psi(\mathbf{r}) = E' \psi(\mathbf{r})$$ where we did the change of variable $(x,y,z) \rightarrow (a_0 x, a_0 y, a_0 z)$, with $a_0 \equiv 1$ Bohr, $E'= ...


12

The first thing to note is that the correspondence between finding roots of a polynomial (any polynomial) and finding the eigenvalues of an arbitrary matrix is really direct, and it's a rich subject, see Pseudozeros of polynomials and pseudospectra of companion matrices by Toh and Trefethen and the references there. Basically, the 2×2 case is trivial and ...


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

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

Given a nonlinear eigenvalue problem of the form $A(\lambda)x = 0$, reducing it to a real equation $\det(A(\lambda))=0$ is known to be a poor method for just the reason you've discovered yourself. The determinant is too ill-behaved for this to work on non-trivial problems. Even for an ordinary eigenvalue problem $Ax=\lambda x$, if you have a good ...


9

Given $A \in {\bf S}^n$ (a positive definite matrix) with eigenvalues $\lambda_1 \leq \lambda_2 \leq \ldots \leq \lambda_n $, then: $\displaystyle f_k(A)=\sum_{i=1}^{k} \lambda_i$ is concave. Why? $$f_k(A) = \inf \left\{ {\bf tr}(V^T A V) | V \in {\bf R}^{n \times k}, V^T V = I \right\}$$ This follows from the Poincare separation theorem (see e.g. Horn ...


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

The simple answer is that you would use inverse iteration (subspace or with deflation). This is basically the power method (repeatedly multiplying the matrix by a vector and normalizing, singling out the eigenvector corresponding to the eigenvalue of largest magnitude) applied to $A^{-1}$. Since you desire the $k$ closest to the origin, you need to use some ...


8

This is trying to compute the eigenvalues by computing the roots of the characteristic polynomial. In this case, the characteristic polynomial is $p(t) = t^3-2t^2x$, $x=1.25\times 10^6$, and zero is a root of multiplicity two. There is generally no good way of computing double roots of a polynomial to precision better than $O(\sqrt{\epsilon})$, because if ...


8

You can use the shift-invert spectral transform [1] and compute the spectrum band by band. The technique is also explained in my article [2]. Besides the implementation in [1], an implementation is available in C++ in my Graphite software [3] (update Jan 17: now everything is ported to geogram/graphite version 3), that I used to compute the eigenfunctions ...


8

There are two relatively convenient options for calculating selected (e.g. a few largest or smallest) eigenvalues using Eigen. The first is Spectra, a header-only C++ library based on Eigen that uses algorithms similar to ARPACK (implicitly-restarted Arnoldi) to calculate a few eigensolutions. Since it is header-only, you simply download and include the ...


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

Note that, if $D$ is invertible, the eigenvalues of $A$ and $DAD^{-1}$ are the same. You can avoid floating-point underflow when forming the matrix by scaling the companion matrix by a diagonal matrix $D$ such that $D_{ii} = 1/\sqrt{(n-i)!}$. For the polynomial $x^4 + a x^3/\sqrt{1!} + b x^2/\sqrt{2!} + c x/\sqrt{3!} + d/\sqrt{4!}$, this gives a modified ...


7

Conclusions are indeed possible, but Gershgorin's circle theorem must be supplemented with other results. Let \begin{equation} \lambda_1 \leq \lambda_2 \leq \lambda_3 \leq \lambda_4 \end{equation} denote the eigenvalues of $M$. Let $N$ denote the lower 3 by 3 corner of $M$, i.e, \begin{equation} N = \begin{bmatrix} 3 &-1 & -1 \\ -1 & 2 & ...


7

There exist special techniques for updating the eigen-decomposition of time-dependent covariance matrices. Given a "prior" eigenvalue decomposition (say at some initial time $t^0$), these recursive algorithms lower the complexity of the spectrum update from $\mathcal{O}(N^3)$ (essentially the cost of a new eigendecomposition) to $\mathcal{O}(k N^2)$ where $N$...


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


7

Unfortunately, convergence of GMRES does not have a clear dependence on the distribution of eigenvalues. It was proved by Greenbaum, Ptak and Strakos in 1996 that you can construct examples with an arbitrary spectrum and an arbitrary convergence history: that is, give me any $n$ nonzero complex numbers, and any decreasing sequence $\|r_k\|$, and I can ...


6

Reformulating (1) and (2), the system reads $$ \begin{pmatrix} \mathcal{L}&\\ &\mathcal{L} \end{pmatrix} \begin{pmatrix} W\\ h_2 \end{pmatrix} = \kappa \begin{pmatrix} -\mathcal{L} &\mathcal{M}\\ \mathcal{M} & -\mathcal{L} \end{pmatrix} \begin{pmatrix} W\\ h_2 \end{pmatrix} $$ where the calligraphic symbols are composed of some linear ...


6

I was able to solve your example by using one of the eigs options to increase the number of Lanczos vectors used in the eigensolution from the default (20 in your case) to 50: opts.p = 50; e = eigs(A,10,'LA',opts) I don't believe the algorithm is extremely sensitive to this number but, obviously, some experimentation with your real problem will be ...


6

@VictorMay has provided an answer using the inverse power iteration, but this is of course expensive. If you have an estimate that $\bar\lambda > \lambda_i$, i.e., it is an upper bound to all eigenvalues including the positive ones, then $A-\bar\lambda I$ is a negative definite matrix. You can then apply the power iteration (instead of the inverse power ...


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

We have the matrix $A$ that can be expressed as $A = JGJ^T$. The first thing is to calculate the QR decomposition of matrix $J$. Because of the low rank of the matrix it can be done very fast with, for instance, modified Gram Schmidt algorithm. Now we can write $A$ as $A = QR G R^TQ^T$, where $Q$ is an orthonormal matrix ($Q^T Q = I$). We define $F$ as ...


6

One possibility is to use a combination of ARPACK and ViennaCL: ARPACK is an eigensolver. It works with a callback interface (you supply a function that computes $Ax$ for a given $x$ and it computes the eigenvalues by invoking this function multiple times). If you want to solve a generalized eigenvalue problem ($ Ax = \lambda Bx$) you need to provide a ...


6

The QR/Francis algorithm is the go-to choice for dense eigenproblems, but there are a few competitors around: The Jacobi algorithm (like QR, another algorithm with an unfortunate name, which can be confused with a method to solve linear systems...). Main idea: choose a nondiagonal element of $A_k$; apply a suitable 2x2 Givens rotations $A_{k+1}=Q_kA_kQ_k^*$ ...


6

This is sometimes known as Buzano's inequality (http://www.jstor.org/stable/2159168). In general, if $\|x\|=1$, $P=xx^{\top}$ is a projection operator, so a simple application of Cauchy-Schwarz leads to $$ 2|b^{\top}xx^{\top}a|-|b^{\top}a| \leq |b^{\top}(2P-I)a| \leq \|a\|\,\|b\|, $$ which gives Buzano's inequality $$ |b^{\top}xx^{\top}a| \leq \frac{\|a\|\,\|...


6

Seems that you have a duplicate eigenvalue. Thus, you have two eigenpairs $(\lambda_1, x_1)$ and $(\lambda_2, x_2)$ where $\lambda_1 = \lambda_2$. Denote $\lambda = \lambda_1 = \lambda_2$. Let $\alpha$ and $\beta$ be arbitrary complex numbers. Then $$A (\alpha x_1 + \beta x_2) = \alpha A x_1 + \beta A x_2 = \alpha \lambda x_1 + \beta \lambda x_2 = \lambda (\...


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