Hello.

Are there any plans to implement the generalized schur decomposition in eigen? How difficult do you think this is?

And what is the most numerically stable way to get generalized eigenvalues with eigen?

Thank you

For selfadjoint problems there is already a GeneralizedSelfAdjointEigenSolver class. For general matrices, there is no plan, and honestly I haven't looked at how to do it robustly/efficiently.

Nevertheless, to solve for:

A x = l B x

with B invertible, you could also compute Q = B.lu().solve(A);

solve for Q y = l y

and if you need the eigenvectors with appropriate norms:

x = y / (y' B y);

Thnk you!

I am looking for (any) real eigenvalue l of the generalized eigenvalue problem Ax=l Bx, where A and B are real symmetric matrices, but often neither one is positive definite...

that's more difficult if they are both singular.

Well, if any one is singular, then l=0 or l=infinity is my answer... So it's safe to assume that they are non-singular and symmetric, but not necessary positive definite.

OK, so the second option I proposed should still work: compute Q = B.lu().solve(A); and then solve for Q y = l y should do the job. It is not the best, but the simplest.

Thank you.

I found a nice article on how to compute generalized Schur decomposition in Golub, Van Loan, "Matrix Computations". I would like to try implementing it for Eigen, if somebody could guide me.

Sure, we can guide you. You could start to have a look at the RealSchur.h file to see how much code can be exploited for the generalized version. The RealSchur implementation is a bit complicated because it strives to avoid complex arithmetics.

I am writing QZ decomposition as in the book, without trying to optimize anything yet. I have a couple of questions:

1. How do I check if a square middle block of the matrix is singular? Block, as well as the whole matrix, is upper-diagonal. Should I use the condition number of A.block(k,k,m,m)? What is eigen code for that?
2. At some point I have to assume that a(i,i-1) is zero, if |a(i,i-1)| < eps ( |a(i-1,i-1)| + |a(i,i)| ), where eps is a small number... what small number should I use?

Another note: typical step of qz algorithm is: compute Householder reflection Q such that (a, b, c) Q = (0, 0, *), and apply it on the right to 3 columns in a big matrix. Or the same for two-dimensional Q. I guess I will need a pretty efficient implementation of this step. Any advice here?

We have everything you need to construct and apply householder reflectors in the Householder module. Usage examples in HouseholderQR, Tridiagonalization, and HessenbergDecomposiiton.

To detect x can be considered to be zero, you need a reference non zero value ref, and then call:

if (internal::isMuchSmallerThan(x,ref)) ...

The choice of ref depends on the application. For instance, to detect a triangular matrix B is singular you can something like:

if(internal::isMuchSmallerThan(B.diagonal().cwiseAbs().minCoeff(),B.diagonal().cwiseAbs().maxCoeff())) {//singular}

Thank you!

I've got a first working version of QZ decomposition for real numbers. It uses double implicit shift, as in "Matrix Computations". This seems to be a good algorithm, at least I don't know any better.

I don't understand how to use Householder module for my needs, I couldn't find an example that would be clearly stable and fast. I am computing reflections matrices explicitly... I would very much appreciate some clarification here. Typical problem is in my post above.

I spent a lot of time debugging, now going out for a ride.

It seems your usage of Householder reflectors is unusual, since in the three examples I listed, we use n-D reflectors such that:

(a,b) * Q = (*,0)

where a is a scalar and b a vector. b is called the essential part.

In your case you are specifically interested in 3D (and 2D?) reflectors so maybe our general Householder methods are not ideal here. Nevertheless you can do:

Matrix<Scalar,3,1> v(c,b,a);
Scalar tau;
RealScalar beta;
v.makeHouseholderInPlace(tau, beta);

Now: (c,b,a) * H == (0,0,beta)

Matrix<Scalar,Dynamic,1> ws(m*2); // workspace (should be declared once)
mat.block(i,j,m,3).transpose().applyHouseholderOnTheLeft(v.template tail<2>(), tau, ws.data());

Hope that helps.

Thank you again.

Yes, I do "zero chasing", which seems to be the most time-consuming part of procedure. I am applying the Householders reflectors from the right to change a piece of matrix (A(n,m), A(n,m+1), A(n,m+2)) to either (0,0,*) or (*,0,0), or Householder reflectors from the left to do the same on subcolumns
(
A(m,n)
A(m+1,n)
A(m+2,n)
),
or the same for subrows or subcolumns of length 2.

One more question:

I have two arrays to sum up, one with 7 coefficients and the other with 11. To minimize roundoff errors I store them in std::vector, then sort them by absolute value, and sum up from the smallest values. What would be a proper way to do this? My current code is

Code: Select all

bool scalar_abs_comp(Scalar x, Scalar y) {
    return std::abs(x)<std::abs(y);
}
...
                std::vector<Scalar> xx(11), yy(7);
                Scalar x=0,y=0,z;
                xx[0] = a11*a11*b11i*b11i;
                xx[1] = a12*a21*b11i*b22i;
                ... (silimar expressions for all of them) ...
                yy[6] = a21*a98*b89*b11i*b88i*b99i;
                std::sort(xx.begin(), xx.end(), scalar_abs_comp);
                for (Index it=0; it<xx.size(); it++)
                    x+= xx[it];
                std::sort(yy.begin(), yy.end(), scalar_abs_comp);
                for (Index it=0; it<yy.size(); it++)
                    y+= yy[it];