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


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)


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

The trick is trying to find out why that matrix has real eigenvalues in the first place. Usually it is because a suitable set of conjugations turns it into a symmetric matrix, and then you can reduce to a symmetric computation. Multiplying and dividing by $(AC)^{-1}$ you can rewrite $$ D_1 = C^{-1}BAC (C^{-1}BAC+I)^{-1}, $$ so your computation is equivalent ...


5

If we assume that $D$ is nonsingular, then there is a relatively straightforward (and efficient) solution based on an $LU$ decomposition. If we write $$ \pmatrix{D & B \\ B^T & A} = \pmatrix{ L_{11} & \\ L_{21} & L_{22}} \pmatrix{U_{11} & U_{12} \\ & U_{22}} = \pmatrix{L_{11} U_{11} & L_{11} U_{12} \\ L_{21} U_{11} & L_{21} U_{...


5

You could define a linear opearator and pass it to the function eigsh. Ideally, your matrices $L$ and $X$ are sparse so you can take advantage of the matrix-vector product. In your case, you would have something like the following. import numpy as np from scipy.sparse.linalg import LinearOperator, eigsh def mv(v): a = 2.3 return L@(L@v) + a*X@(X.T@v)...


5

The other answers already tell you what went wrong, but I will add a terminology note: the term for what is happening is that the pencil $A - \lambda B$ is a singular matrix pencil, i.e., $\det (A - \lambda B)$ is identically equal to zero. So there are no generalized eigenvalues (or, at least, they cannot be defined as usual as the roots of the generalized ...


4

Looks like everything is working relatively OK? Your matrix is of order 1e10, so residuals of 1e-4 are actually close to machine precision. The convergence criterion is indeed violated, but not by much; not sure what's going on there, are you sure Arpack really guarantees it or is it a best effort kind of thing? It surprises me that you get different results ...


3

Yes, the eigenvectors found with this method may depend on $x$ and $y$, but no, it doesn't matter in practice. If $A$ and $B$ share a basis of common eigenvectors, then $$ A = V\operatorname{diag}(\alpha_1,\dots,\alpha_n)V^{-1}, \quad B = V\operatorname{diag}(\beta_1,\dots,\beta_n)V^{-1}. $$ If $A$ has all distinct eigenvalues, then $V$ is unique, up to ...


3

I suspect the root of your trouble is what has been detected in the comments by Vibe: For any number $\omega\in \mathbb{K}$ (with $\mathbb{K}= \mathbb{R}$ or $\mathbb{C}$) you can find $\boldsymbol{X}$ such that $AX = \omega BX$ (with $A$ and $B$ taken in your concrete example). You have already decomposed the problem in 4 blocks of 3 variables. Then let us ...


3

Using the Cholesky decomposition is the quickest way I know of checking if a symmetric matrix has negative eigenvalues. Nothing wrong with that. Plus, if it succeeds, you already have the Cholesky decomposition! Of course, if there are any negative eigenvalues, it will fail. You will not need to check the entries of $L$.


3

Cholesky factorization is for symmetric positive definite matrices, and it will fail if the matrix has negative eigenvalues. You should use singular value decomposition for that purpose, or maybe a QR algorithm would suffice if you just need some of the eigenvalues. Edit 1: Of course, I should also add that, in general, L would not give you much information ...


2

It's a bit late, but I have a very simple answer: there is nothing wrong with the code I was simply missing the extra optimisation that Newman recommends in his paper. The results are perfectly in-line with what is reported in the paper without applying the extra optimisation step. I'll leave this up in case it helps someone implementing the algorithm.


2

Your problem was the lower limit of integration. It should have been $-x_e$ instead of 0, since $x_e$ is the equilibrium point for the potential and not the minimum distance. After correcting that, you get the following #%% Solution xe, lam = 1.0, 6.0 # parameters for potential xmax = 10 # Bval, Bval2 = wavefunction values at x = bound1, bound2 bound1, ...


1

OK, so with new version in edit: it's not your fault, it's https://github.com/JuliaLinearAlgebra/Arpack.jl/issues/87. You can either call Arpack manually yourself, or even better use a pure-julia solver (KrylovKit.jl and IterativeSolvers.jl are good choices)


1

Probably you accumulated rotations wrong. At each step, you have a working matrix $D$ that starts from $D=A$ and should converge to a diagonal matrix, and a matrix $Q$ that accumulates all rotations performed on $D$. Print after each step $\|QDQ^T - A\|$, and see where it stops being small.


1

First some comments on why such matrices are hard for solvers. Notice that if $U$ is upper diagonal like [[. 2 . . .] [. . 2 . .] [. . . 2 .] [. . . . 2] [. . . . .]] then $(I - U)^{-1} = I + U + U^2 + ... U^{n-1} \ (U^n = 0)$. For this example, $(I - U)^{-1} =$ [[ 1 2 4 8 16] [ . 1 2 4 8] [ . . 1 2 4] [ . . . 1 2] [ . . . . 1]] ...


1

This answer is a bit late, but I was trying to implement the code you wrote and I found some issues. You have not implemented the Lanczos algorithm correctly. The diagonalized matrix, $T$, itself is correct. So when you try and find the eigenvalues you obtain from the diagonalized matrix they are correct, but the matrix $V$ is not correct.The matrix $V$ is ...


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