Hello everyone!
I am currently working on a generalized eigenvalue problem, defined by the equation Aϕ=λBϕ, where matrix B is singular. Both matrices A and B are defined as sparse matrices, as I aim to apply this script to high-order linear systems.
Matrix B is defined as follows: B = spdiagm(n, n, 0 => ones(ndif)) # with n >> ndif
Can anybody help me calculate the eigenvalues and eigenvectors of this system? I have checked the LinearAlgebra, ArnoldiMethod, and Arpack packages, but I don’t quite figure out the best way to solve this problem.
Thank you very much for your help!
If A is not singular you can move \lambda to the other side. If they are both singular then you may have a “singular matrix pencil” and things get tricky (google it).
2 Likes
Hello!
Thank you very much for your answer; it’s been very helpful. Unfortunately, matrix A is also singular (I am solving a DAE system). Do you have any suggestions for handling this kind of problem?
If both A and B are singular, you might still be in luck. A pencil A-\lambda B is singular if A-\lambda B is singular for every \lambda. (Otherwise it is called “regular”). If there exists \sigma for which A-\sigma B is nonsingular then you can solve the shifted problem
B x = \mu (A-\sigma B) x
for \mu = 1/(\lambda - \sigma) and then you have \lambda = \sigma + 1/\mu
If no \sigma exists that makes A-\sigma B nonsingular your pencil is singular and you have a much, much harder problem. The only way to even meaningfully define eigenvalues at that point requires thinking about the Kronecker canonical form and reducing to a smaller regular pencil. The algorithms for that are direct (not iterative, at least as far as I know), complicated, involve delicate rank decisions, and are not guaranteed to identify Kronecker structure you might be interested in. There are a bunch of references at:Templates. Look at the chapter on singular pencils. I’d explore other options first.
Edit: The comparison I made to computing the Jordan form is perhaps overstating the difficulty since you don’t have to fully reduce to Kronecker form. But it’s still very challenging.
7 Likes
Julia code (or wrappers) for reduction to Kronecker-like staircase form and eigenvalue computation for pencils resembling yours (as dense matrices) is included in MatrixPencils.jl. The specific form of B you mention arises in control theory, so you may also be interested in DescriptorSystems.jl
4 Likes
Thank you very much for your answer! It has been very helpful, but I have encountered a limitation. The system I am representing is a Differential-Algebraic Equation (DAE). This solution works great when the number of algebraic equations is comparable to the number of differential equations. However, when n_dif≪n_alg, I am unable to “regularize” the problem with the shifted problem approach. Since I am working with power systems, I believe this is a common configuration (n_dif≪n_alg). Do you have any suggestions for resolving this issue? Thank you!
Is there a reason you aren’t using a symbolic interface like MTK to dimplify away some of the algebraic equations?
Hello, thank you for your response! Indeed, it is a power systems problem, and I have checked DescriptorSystems.jl. However, as far as I understand, it does not work with sparse matrices. When I try to reduce the system myself, I get blocked when attempting to execute a J1 \ J2 operation with sparse rectangular matrices. If you have any ideas on how to solve this, I would be grateful. 
I’m not sure I fully understand the question, but I believe the answer is that I have the A matrix of the system, which is relatively large. Therefore, I’m unsure if I can use a symbolic interface. To be honest, I haven’t checked this solution, so I will look into it. Thank you!
AFAIK DescriptorSystems should be compatible with sparse matrices. If there is a specific error, you should put it here along with the version of the package you are using and your environment. Or you could open an issue on the github.
1 Like
What do you mean by relatively large ?
1 Like
I’m not sure I understand your problem well enough to offer anything useful. You originally asked about a generalized eigenvalue problem of unspecified size with B singular. Then A was singular as well. I don’t actually understand why the shift and invert didn’t work for your problem, maybe because I don’t understand exactly what you need to do. Or maybe you do in fact have a singular pencil? (That word “regularize” confuses me in this context.) I also am not sure exactly what you want to do with your descriptor system. A more complete statement of the problem might help.
1 Like
As @mstewart said, all the Kronecker-like algorithms typically generate full matrices, so they aren’t formulated for sparse arguments. But it sounds like you don’t need that machinery.
The pencils corresponding to DAEs describing power system dynamics should not be singular. (Otherwise one does not have a consistent DAE - solutions may be nonunique or nonexistent. That could happen if you try to model a disconnected network. This differs from control theory applications which may be rank-deficient and therefore singular by design.) That means the shift-and-invert method @mstewart suggested should work, i.e. that the shifts can be chosen so that the factorization of A - \sigma B is not badly conditioned, even if B has modest rank.
If you can’t get a well-conditioned factorization with several shifts, I can think of a couple possible causes.
You may have (effectively) redundant constraints. This would make a singular pencil, but because of a mistaken formulation rather than physics. From your description, the constraints are in the bottom of A, so perhaps you could do a sparse QR of that part of A to identify the numerical nullspace. Then either swap in new constraints to handle it or drop some spurious degrees of freedom by projecting out the nullspace.
Alternatively, there may be unfortunate scale choices. You may have physical units for different elements of the representation which make it numerically singular when the problem is actually benign. You may be able to diagnose that by trying to integrate the DAE.
3 Likes
I think this is a good point to make, although I didn’t know enough about power systems to try to make it. It’s probably not a very accessible reference outside of buying the book or going to a university library, but if anyone is curious about this, the potential nonexistence and nonuniqueness of solutions to DAEs associated with singular pencils is described in detail in Matrix Theory vol. II by Gantmacher using the Kronecker canonical form.
2 Likes
Hello,
You are right; I will try to define my problem more clearly (apologies for any misunderstanding – I am also discovering more about the problem and its implications as we discuss it!).
- I have a linearized DAE system defined as : B \dot x = A x. I know A and B matrices. I believe this can also be referred to as an implicit linear system or a descriptor system.
- I want to calculate the eigenvalues of this system, which means I need to solve the generalized eigenvalue problem A \phi = \lambda B \phi.
- Since the system involves approximately 20e3 differential variables, I need to work with sparse matrices.
- When I try to solve the shifted problem, it works if rank(A) \approx rank(B), and it fails when rank(A) >> rank(B).
As @Ralph_Smith suggested, I am going to check if my problem is correctly formulated, so that I can solve the shifted problem even if B has a lower rank.
Thank you again for your help; it has been incredibly useful!
Ok, thank you for the suggestion and the clear and complete explanation! I’ll check it out.
I think you’re probably on the right track in following the suggestion from @Ralph_Smith. A singular pencil seems likely to indicate something is wrong with your formulation of the problem. But I did want to add one thing: 20,000 variables is not something that is out of reach for dense eigenvalue solvers if you have a decent amount of memory and some patience. The QZ algorithm (eigen(A,B) in Julia) would be fairly time consuming, but if you get shift and invert working, then the QR algorithm applied to the shift and invert isn’t too terrible. I just did a 20000\times 20000 problem in about 20 minutes on an unimpressive computer. This might be an option if you don’t want to worry about the Arnoldi method or are interested in eigenvalues that are hard to get using using the Arnoldi method. (i.e. those that are in the interior of the spectrum). Although, depending on which shifts result in a reasonably well conditioned matrix, you might be able to shift and invert to make your eigenvalues of interest large so that Arnoldi will find them easily. In any event, if Arnoldi poses challenges, dense methods might be an option.
3 Likes