for such small sizes you should not bother and either directly accumulate your entries, or use fixed size objects:
Matrix<Scalar,7,1> xx;
Matrix<Scalar,11,1> yy;

xx.sum();
yy.sum();

For fixed sizes, Eigen uses a recursive divide and conquer strategie that naturally reduces the risk of precision loss.

btw, do you plane to share your code with us so that we can integrate into Eigen itself? Would be great for Eigen! If so, may I suggest you to create a fork on bitbucket (https://bitbucket.org/eigen/eigen) so that we help you in this task and more rapidly converge to a stable implementation.

Thanks.

Yes, definitely. I'd be glad to see my implementation become part of Eigen. I didn't post any code yet because I am still changing things a lot.

Ok, I will try making fork.

excellent,
I'm looking forward.

ggael wrote:for such small sizes you should not bother and either directly accumulate your entries, or use fixed size objects:
Matrix<Scalar,7,1> xx;
Matrix<Scalar,11,1> yy;

xx.sum();
yy.sum();

For fixed sizes, Eigen uses a recursive divide and conquer strategie that naturally reduces the risk of precision loss.

Just tried summing the values without sorting. On my test case it gives error of about 1e-15, while with sorting only 1e-18 with Scalar=long double.

The problem there is to compute first column of matrix (M-lambda1 I)(M-lambda2 I), where M=A B^{-1}, A is upper Hessenberg, and B is nonsingular upper triangular, and lambda1 and lambda2 are eigenvalues of M's last 2x2 block. That column has only first three non-zeroes x,y,z, and they be computed in O(1). I have derived an analytic solution, but it is quite a large formula. I wonder if there's any better way to do this.

I have forked eigen, see
https://bitbucket.org/Khumarahn/eigen-qz/

I implemented one class RealQZ, and didn't touch anything else. My implementation is quite raw, and prints debug messages.

Where do I put "#include <vector>" and "include <algorithm>" ?

Should I do a check that input matrices are not complex? How?

Suppose I do a Householder QR of matrix B: B = Q R. What is the best way to compute Q* A for some other matrix A?

Usage example:

Code: Select all

        // make two random matrices
        Matrix<long double, Dynamic, Dynamic> A(dim,dim), B(dim,dim);
        A.setRandom(); B.setRandom();
       
        // perform real QZ; A and B are not changed
        RealQZ<m_type> qz(A,B);

        // take matrices
        Matrix<long double, Dynamic, Dynamic> a = qz.matrixA(), b = qz.matrixB(),
           q = qz.matrixQ(), z = qz.matrixZ();
        // now A = q a z, B = q b z,
        // b is upper triangular and a is upper quazi-triangular,
        // q and z are orthogonal


Hi,

thanks for creating the fork. Regarding (M-lambda1 I)(M-lambda2 I), here is how I would compute its first column (only the first 3 coeffs are nonzeros):

[CODE]
Vector2 y(A(0,0)/B(0,0)-lamba2, A(1,0)/B(0,0));
Vector3 z;
z << A.template topLeftCorner<2,2>() * (B.template topLeftCorner<2,2>().template triangularSolver<Upper>().solve(y)) - lambda1 * y, A(2,1)*y(1);
[CODE]

If I did not make mistake, that's not so complicated.


How have you computed the errors? Are they the relative errors of the resulting factorization?

You can compute Q * A as follow:

HouseholderQR<MatType> qr(M);
QA = qr.householderQ() * A;

I'll look at the code later.

Hi.

I am trying to figure out the source of numerical instability. This is how I compute errors: let A and B be two matrices, A = q a z and B = q b z be their qz decomposition. Then I compute Frobenius norms |q* q - I| , |z* z - I|, |q a z - A|/|A|, |q b z - B|/|B|. Here |x| is for Frobenius norm, indeed. May be, I should be computing errors differently?

I work with long double, where epsilon is about 1e-19. I fixed dimension = 5 and run decomposition for random matrices in a loop. All the norms above max at about 5e-18 except |q a z - A|/|A| that reaches 1e-15 sometimes; and (very rarely) for some pairs of matrices the procedure does not converge.

I am about to start investigation.

Strange... At some point I have an unreduced upper Hessenberg A, and I am looking for zeros on its sub-diagonal - if I find any, it is possible to "deflate" - reduce the problem to a lower dimension. In my test example I meet a block on the diagonal:

Code: Select all

      1.7    -0.101
  1.6e-15       1.8


1.6e-15 belongs to a sub-diagonal. If it is "zero", I can deflate. To check this, I call a function

Code: Select all

internal::isMuchSmallerThan(1.6e-15, 1.7+1.8)

It is true. So I set 1.6e-15 to zero. Then this is a biggest part of resulting error of calculation. I replaced check with

Code: Select all

 
internal::abs(1.6e-15) <= NumTraits<Scalar>::epsilon() *(internal::abs(1.7)+internal::abs(1.8))


and accuracy is better now.

For a test example when procedure didn't converge, I observe a cyclic behavior. It probably can be taken care of with some random shifts.

UPD: With all the last modifications, including random shifts on iterations 15, 31, 47, etc., I have solved more than 300000000 qz problems by now, and all converged, and maximal error was about 6.7e-18, which looks acceptable. Algorithm makes about 9.5 iterations on average... I might think about some optimizations and code style now...

UPD: With all the last modifications, including random shifts on iterations 15, 31, 47, etc., I have solved more than 300000000 qz problems by now, and all converged, and maximal error was about 6.7e-18, which looks acceptable. Algorithm makes about 9.5 iterations on average... I might think about some optimizations and code style now...

great work!!

I updated the program. I almost completely rewrote it, and it looks much like what I found in RealSchur. Take a look at my fork
https://bitbucket.org/Khumarahn/eigen-qz/

I have found some cases where program would not converge, and I am still investigating this, and then there will be some minor things to change.

Would be good if somebody took a look at my implementation and let me know if something has to be changed. May be we should continue discussion in mailing list instead of this forum?

I cleaned up, did what was left, and am nearly done.

Should I use A.coeffRef(i,j) instead of A(i,j)? Or A.coeff(i,j)? Is there a real difference? And I probably should replace everything like

Code: Select all

0.5*(a+b)

by

Code: Select all

Scalar(0.5)*(a+b)

Hi, sounds pretty good. Some comments though:
- yes to add explicit conversion to Scalar(.)
- coeff() / coeffRef() skip the checking of the indices wrt matrix sizes, so they are recommended for Eigen's code.
- I did not looked closely at the algorithm, but it seems that in splitOffTwoRows, givens rotations that are applied to S and T could be optimized by taking into account that they already contain many zeros, isn't it? Same in pushDownZero?
- The most important for an immediate inclusion is that it lacks a unit test.
- Many block operations can be optimized, e.g.:
- m_T.block(0,k+2,lr,1) -> M_T.col(k+2).head(lr)
- m_T.block(0,k,lr,2) -> Block<MatrixType,Dynamic,2>(m_T,0,k,lr,2) or m_T.middleCols<2>(k).topRows(lr)
- some for m_S.block(k,fc,3,dim-fc)
- minor: we usually use only 2 spaces for the indentation

ok, thank you for comments

What is the difference between A.coeff(i,j) and A.coeffRef(i,j)? This is probably explained in documentation, but I couldn't find it.

What is the unit test? Results of testing the algorithm on random matrices?

Yesterday I realized there are algorithms that do QZ decomposition of large matrices much faster than my program...

coeff() returns by value or const ref when possible. coeffRef returns a writeable reference.

You can have a look at our unit test in test/. For instance, see test/eigensolver_generic.cpp. If our file is named qz_real.cpp, then it must implement a
void test_qz_real();
function that does the tests on random matrices of different sizes and scalar types.

Regarding performance, would it be possible to exploit our HessenbergDecomposition algorithm instead of doing it via Givens rotations?