QR update

[peng.du]

Hi there

Is the matlab qrupdate functionality implemented in LAPACK or anyone working on it?

Thanks.

[Julien Langou]

Hi Peng, as far as I know nobody in the team is working on or planning to work on updating the QR factorization of a matrix after a rank-1 update of the coefficient matrix (Mathworks' "qrupdate"). The functionality was in LINPACK and never made it to LAPACK. Contributions welcome in that respect though ;) Cheers, Julien.

[sven]

Daniel Kressner has some code at:

www.math.ethz.ch/~kressner/qrupdate.php

Hopefully there will be some LAPACK code in the future.

Best wishes,

Sven Hammarling.

[peng.du]

thanks guys. I'll see what I can do.

BTW, awesome report: http://eprints.ma.man.ac.uk/1192/01/cov ... 08_111.pdf

[peng.du]

So the idea is more or less the same: using givens to zero out spikes. The original Q is taken from the result of DGEQRF in the form of elementary reflectors, and at a certain point the final Q is generated by applying a series of givens matrices to Q from left or right. This makes the new Q in its full form. Sorry for being lazy but any suggestion on routines that squeeze it back into the elementary reflectors form?

[sven]

Daniel's routines only use elementary reflectors.

Sven.

[bravegag]

Hello,

I have the same problem as peng.du has, trying to reuse Daniel's routines I pass as input as A the old QR as returned by DGEQRF and B is a matrix that contains the extra row and doesn't work for me.

My use-case is simple I create a QR factorization using DGEQRF, then I have a pointer to the QR compact LAPACK form which I can pass around.

Daniel's A input is a squared upper triangular matrix ... I am confused by A, can anyone clarify?
The B matrix input I think I understood from his test case to be a matrix with the extra block of rows to be updated.

My C code look like this trying to append 1 row (I renamed his dbqru.f to addrows.f):

extern void F_addrows(lapack_int* m, lapack_int* n, double* A, lapack_int* lda, lapack_int* b,
lapack_int* ldb, double* tau, double* work, lapack_int* lwork, lapack_int* info);

lapack_int m = 1;
lapack_int n = B->rows - 1;
lapack_int lda = n;
lapack_int ldb = m;
lapack_int info = 0;

// use as temporary memory
int lwork = n*n;
tvector* work = vector_init(lwork);
vector_realloc(B->tau, MIN(B->rows, B->cols));
F_addrows(&m, &n, B->QR->data, &lda, B_extra->data, &ldb, B->tau->data, work->data, &lwork, &info);
if (info != 0) {
// QR update failed
fprintf(stderr, "F_addrows failed with info = %d\n", info);
exit(EXIT_FAILURE);
}
vector_delete(work);

Can anyone advice?

TIA,
Best regards,
Giovanni

[bravegag]

Hello,

I have been trying to reuse Daniel's routine:
http://www.math.ethz.ch/~kressner/qrupdate/dbqru.f

Can anyone please advice how to update QR from Daniel's output? given that the old QR is in compact LAPACK form as result of calling DGEQRF, meaning the upper triangle contains R and the lower triangle a normalized set of reflectors and tau. I can call Daniel's routine fine and I get the B replaced with column reflectors Hnew_1 ... Hnew_n ... and get the tau_new. My problem is how to _merge_ these new smaller reflectors with the existing ones Hold_1 ... Hold_n and tau_old in a way I can put them all as one back in compact form as if it was output by DGEQRF.

I know I can do this by generating Q explicitly first and then compacting it all back ... but it is a lot of unnecessary flops while I only need to merge the existing and new reflectors.

Many thanks in advance,
Best regards,
Giovanni

[Julien Langou]

Hi Giovanni,

First of, it is great that you are using Daniel's subroutines but maybe you can try to use the LAPACK 3.4.0 ones. The subroutine name is DTPQRT.

So consider A1 of size M1-by-N with leading dimension LDA1 and A2 of size M2-by-N with leading dimension LDA2.
We also assume that

Code: Select all

M1 >= N

(So as to be able to store a full N-by-N triangle in A1.)

First perform the QR factorization of A1 with DGEQRF. (You may also use DGEQRT or DGEQRT3 actually.)
Then perform the updating. So

Code: Select all

DGEQRF( M1, N, A1, LDA1 ... )
DTPQRF( M2, N, 0, NB, A1, LDA1, A2, LDA2, ... )


Note: NB is the blocking size of the algorithm. We decided to have it part of the interface so you (the user) need to decide what you want NB to be.
In LAPACK, in general, NB is set by the ILAENV( ) function. So this kind of interface slightly differs from what is standard in LAPACK.
We have NB in the interface of DGEQRT and DTPQRT (both introduced in LAPACK 3.4). We also decided to export the " T matrix " to the user, see the
T, LDT arguments in the interface. This is useful for performing efficiently the updates. (This avoids to have to recompute T.)

At this point, you can (as you did with Daniel's code) check that the R factor (stored in A1 as the upper N-by-N part) matches the one of DGEQRF with [ A1 ; A2 ]
up to signs in the rows of R.

Note you can also do:

Code: Select all

DGEQRF( M1, N, A1, LDA1 ... )
DGEQRF( M2, N, A2, LDA2 ... )
DTPQRF( N, N, N, NB, A1, LDA1, A2, LDA2, ... )


So QR factorization of A1, QR factorization of A2 and then reduce the two triangles to one.
You would get more parallelization out of this. (Less locality.)
This is the same number of operations overall.

You can also do a lots of other things. I let you imagine.

Can anyone please advice how to update QR from Daniel's output? given that the old QR is in compact LAPACK form as result of calling DGEQRF, meaning the upper triangle contains R and the lower triangle a normalized set of reflectors and tau. I can call Daniel's routine fine and I get the B replaced with column reflectors Hnew_1 ... Hnew_n ... and get the tau_new.

( All your wrote makes sense. )

My problem is how to _merge_ these new smaller reflectors with the existing ones Hold_1 ... Hold_n and tau_old in a way I can put them all as one back in compact form as if it was output by DGEQRF.

The answer is that with Daniel's subroutines or LAPACK-3.4, you cannot compute back the exact same reflectors as DGEQRF would do.
The main point is that there is not much point in computing these reflectors. What you want to do is being able to apply them to another
matrix (and this is done easily by DORMQRT and DTPMQRT) or you want to compute Q (and this is done as well by applying the reflectors
with DTPMQRT to the identity matrix.)

I hope this makes sense. If not we can give a sample example.

Cheers,
Julien.

[bravegag]

Hi Julien,

Thank you so much for your answer it will take me a bit to digest it though because I am new to BLAS and LAPACK and also to Fortran, so is a lot to catch up with. If you could provide examples that would be most welcome and greatly appreciated I am very curious how it would compare performance-wise to Prof. Kressner's implementation.

Meantime I understood that merging the old and new reflectors is not possible, like it is not possible (and not efficient) to compute the explicit new Q and try to compact it as a lower triangular set of new reflectors. My target for pursuing the compact form was to 1) keep code simple by managing only one type of representation [LAPACK compatible], the algorithm I am trying to implement is already quite complex 2) avoid bookkeeping a potentially large and unbounded number of reflectors. But for the sake of efficiency I think the choice of bookkeeping all the small reflector updates is the most efficient solution. I then need to call many times Prof. Kressner's DBQRU_APPQ that applies the small set of reflectors in compact form corresponding to one update step to a given matrix before solving the system something like H1_u3*H2_u3*...*Hn_u3*H1_u2*H2_u2*...*Hn_u2*H1_u1*H2_u1...*Hn_u1*H1*H2*...*Hn*c where ui is the QR update number.

The algorithm I am working on is iterative and keeps mutating the main matrix in multiple ways: adding rows, adding cols, deleting cols and multiple solves in-between. So the efficient use of QR updates is critical.

Many thanks for your help.
Best regards,
Giovanni

PS: Are you one of the authors of this paper "Parallel tiled QR factorization for multicore architectures"? I had to do a seminar about it in my University ... it was a really hard cookie to understand but very interesting

[bravegag]

Hi,

And another question, were the addcols and delcols functions also integrated into LAPACK 3.4.0?
http://www.maths.manchester.ac.uk/~clucas/updating/

These functions AFAIK were implemented by the authors of "Updating the QR factorization and the least squares problem, with Sven Hammarling, MIMS EPrint 2008.111. November, 2008."

If not, is there any LAPACK 3.4.0 alternative to the addcols and delcols working?

Many thanks in advance,
Best regards,
Giovanni