Hot answers tagged

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


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


6

LAPACK has been on the cutting edge for just about three decades, and probably still is for its niche. However, given given recent developments in libraries for the simpler BLAS-type matrix operations that LAPACK traditionally builds upon, it is perhaps conceivable that we could see the emergence of serious competitors to the traditional FORTRAN-based LAPACK ...


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

And just a few minutes after asking I found an answer. The procedure above is called "updating LU". This question has a nice generic answer with links to other more specific questions. full rank update to cholesky decomposition So the short answer is no, in my case. You cannot update LU decomposition with a full rank matrix in less than $O(n^3)$, ...


4

As far as I know, there are no such methods in LAPACK. Since LAPACK is the linear algebra package, no nonlinear solvers are included. However, you can use the underlying BLAS for implementing iterative methods. For nonlinear solvers using Jacobians (e.g., Newton's method), the matrix factorizations of LAPACK may come in handy. You may want to have a look at ...


4

LAPACK doesn't have a specialized routine for computing the eigenvalues of a unitary matrix, so you'd have to use a general-purpose eigenvalue routine for complex non-hermitian matrices. This is slower than using a routine for the eigenvalues of a complex hermitian matrix, although I'm surprised that you're seeing a factor of 20 difference in run times. ...


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

The reference LAPACK library prioritizes portability over performance, so it's missing many optimizations, such as vectorization and threading. (Additionally, it relies on BLAS for much of it's performance, so using the reference BLAS will significantly limit performance.) However, most LAPACK users don't use the reference library, but various optimized ...


3

There is another factorization you could consider: the Hessenberg upper-triangular reduction. It's usually used as a preprocessing step in the QZ algorithm, but it has other uses as well. Consider the reduction $(A,M) = Q^T(H,T)Z$. Then $\frac{M}{\mu_j} + A = Q^T(\frac{T}{\mu_j} + H)Z$, which is relatively cheap to solve as it only involves orthogonal ...


3

LAPACK doesn't include any iterative solvers. The routines in LAPACK are for eigenvalues, matrix factorizations, and solutions of systems of equations involving dense matrices while iterative methods are generally used for matrices that are large and sparse.


3

It would help if your question included the matrix written down in the usual format. In addition, it looks like there is a typo in your equation (the last $\sigma_i$ should be $\sigma_{i+1}$?) Anyhow, if I guess correctly what you mean, your matrix is weakly diagonally dominant and irreducible, so by a corollary of Gershgorin's circle theorem 0 is an ...


2

Fortunately, LAPACK provides routines to deal with the $\mathbf Q$ factor from the $\mathbf A = \mathbf Q \mathbf R$ decomposition, [dgeqrf]. To find the projection of an arbitrary $\mathbf B$ onto the space orthogonal to $ \mathrm {range}(\mathbf A)$, you want to form $\mathbf C = \left(\mathbf I - \mathbf Q \mathbf Q^T\right) \mathbf B$. Here are two ...


2

Base case: Computing the factorization requires $\frac 2 3n^3$ operations and the inverse requires $\frac 4 3n^3$ operations. Computing the trace adds $O(n)$. Let's round that to $2n^3$. LU case: Computing the factorization requires $\frac 2 3n^3$ operations. If you were to compute $L^{-1}$ and $U^{-1}$, those operations each require $\frac 1 3n^3$ ...


2

Another approach, which might be of interest to you is randomized sampling. This is of particular interest if you can quickly compute matrix-vector products $x\rightarrow Ax$ and $x\rightarrow A^* x$. The core idea is to form a small sampling matrix $S = A\Omega$, where $\Omega$ is a Gaussian random matrix. If the sampling matrix is large enough, $S$ will ...


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