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

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

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

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

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


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


4

I don't think there is a way to display the used LAPACK function names natively during runtime using the interfaces provided by scipy.linalg. Depending on your goals you can: read the source code and deduce the logic from there. Unfortunately, this is not a runtime-use scenario, but a human analysis. fork your own version of scipy, add custom outputs (to ...


4

You might want to use the fact that: $$ ||A||_2=\sigma_\max(A) $$ where $\sigma_\max$ is the largest singular value. If you are interested in details, this Math SO question should be interesting. Thus, $$ ||A^{-1}||_2=\frac{1}{\sigma_\min(A)} $$ where $\sigma_\min$ is the smallest singular value. You certainly want to avoid the actual calculation of the ...


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

Test your code. "I don't know if my linear algebra routines work or not" is a problem you can solve easily. It is easy to check if you have computed the correct eigenvalues or not; just check that $VDV^{-1}=H$, or if you don't fancy the inversion $HV=VD$. If you are not sure that your individual pieces work, you are walking in the dark. Possibly a hot take: ...


2

If you're lucky enough to have a complex-hermitian $\mathbf A$, the eigendecomposition of $\mathbf A$ can be computed using basically the same algorithmic machinery as the real-symmetric case: an initial/frontend pass to reduce to tridiagonal form $\mathbf T$ via Householder reflections, followed by shifted-QR iterations on $\mathbf T$. Notably, a complex-...


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

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

Your whole code is not understandable to me. Try this alternative components: int grid(int i, double a, double h) { return a+i*h; } int index(int i, int j, int N) { return i*N+j; } //potential function. returns the integral over the grid square indexed by i double V(int i, double x0, double hx, int j, double y0, double hy) { double xi=grid(i,...


1

I believe it can be done with a semidefinite program by adding a multiplicative slack variable. Basically, $$ \begin{array}{rcl} \min\limits_{A \in \mathbb{R}^{n \times n}, P\in \mathbb{R}^{n\times n}} &&f(A)\\ \text{st} && b I \preceq PA + A^TP \preceq a I\\ && P \succ 0 \end{array} $$ Essentially, this is the Lyapunov stability ...


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


1

It turns out the issue was very simple (and not related to a limitation of numpy's matrix power function as such). I initially thought there was some numerical floating point error being propagated - but the example matrix I was testing on contains only integers. The problem was that at $k=20$ the value of the matrix entries exceeded numpy's maximum possible ...


1

The Frobenius norm is not an operator norm, it is a norm on the vector space of linear operators/matrices, which is not the same thing. Just change it to any other preset norm and it should work. It is also the case that your method of computing matrix powers is not stable. The algorithm used in Numpy is basic repeated squaring, which has no normalization ...


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