# Tag Info

23

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

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

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

13

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'= ... 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 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 ... 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 ... 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 ... 11 You have an ill-conditioned eigenvalue problem. Consider a perturbation$\delta A$in your matrix$A$—with double-precision floats this is$O(10^{-16})$. It turns the eigendecomposition from $$X^{-1}A X = \Lambda$$ into (approximately keeping$X$fixed) $$X^{-1}(A+\delta A)X = \Lambda + \delta \Lambda.$$ This means that the change in eigenvalues is ... 10 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 ... 9 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 ... 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.... 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 ... 9 Unfortunately, I don't think there is a good algorithm to do this efficiently. Given the eigendecomposition$\mathbf A = \mathbf X \mathbf D \mathbf X^T$, one is tempted to project$\mathbf v$onto the eigenvectors by introducing the vector$\mathbf u = \mathbf X^T \mathbf v$, forming$\mathbf A + \mathbf v \mathbf v^T = \mathbf X \left(\mathbf D + \mathbf u ...

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

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

"Get more RAM" may be one of your best options. :) Prices are reasonably low right now, and it's one of the best upgrades you can gift your computer anyway. 10k x 10k is borderline but still doable on modern computers: that matrix takes $10^4 \times 10^4 \times 8$ bytes, that is, 760 MiB. On my laptop that code runs without problems. Another option is ...

8

I think the method has too much implementation complexity and too narrow applicability to be worth it. Though the paper is correct to point out the importance of solving the tridiagonal-symmetric eigenproblem in the course of solving the general-symmetric eigenproblem, it neglects to mention that the "frontend" procedure between these two scenarios (...

8

At least for the second question the answer is yes. See for example Mattheij, Robert MM, and Gustaf Söderlind. "On inhomogeneous eigenvalue problems. I." Linear Algebra and its Applications 88 (1987): 507-531, page 516. (The optimality conditions of your problem, $$Ax+b+2\lambda x = 0 \\ x^Tx = 1$$ constitute an inhomogenous eigenvalue problem)

8

The matrix B (M in the documentation) needs to positive definite according to the documentation: "If sigma is None, M is positive definite", this is in addition to the first requirement "M must represent a real, symmetric matrix if A is real" which your B follows. The eigenvalues of your current matrix B are -1, 1 and 6. So matrix B is ...

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 $$\lambda_1 \leq \lambda_2 \leq \lambda_3 \leq \lambda_4$$ denote the eigenvalues of $M$. Let $N$ denote the lower 3 by 3 corner of $M$, i.e, 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

When one says an algorithm is of order $O(n)$, that may mean that the complexity is given by: $c + b*n$. With every new element you add you increase in runtime (effectively). What mathematically minded people often forget is that these statements do not include how large the constants are. That of course carries over to $O(n²)$ and such. I can not answer ...

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

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

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

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