help-octave
[Top][All Lists]

## Re: matrix_type problems

 From: Vic Norton Subject: Re: matrix_type problems Date: Sun, 3 Dec 2006 17:57:42 -0500

```On 12/3/06, at 3:29 PM +0100, David Bateman wrote:
> As for the test that a matrix is upper triangular, A =
> matrix_type(A,'Upper') will force the matrix A to be of type 'Upper',
> without having to have the type probed. This gives a slight gain in
> performance and is valid for the output of such things as a Cholesky
> factorization. However, this isn't want you wanted. To probe whether a
> matrix is upper triangular something like
> strcmp(matrix_stype(A),'Upper') would be more appropriate.

In fact, David, I do want to force the matrix to be of type 'Upper' if that
will speed up performance.

Here is my situation. I have an iterative routine that updates QR
factorizations of submatrices of an m-by-n matrix A. At each iteration I have
an index vector J of length nJ <= min( m, n ), and I need to compute
e = A( : , J ) \ u;
This can be done with 2 * nJ * nJ * ( m - nJ/3 ) flops. (flop = floating point
operation.) The vector u is fixed, but J and nJ change from one iteration to
the next, usually by the addition or deletion of a single index. The matrix A(
: , nJ ) always has rank nJ.

In my algorithm I maintain a QR factorization of A( : , J ) in two arrays, QQ
and RR, of sufficient size. Thus A( : , J ) = QQ( : , 1:nJ ) * RR( 1:nJ , 1:nJ
), so that e can be computed as
uQ = QQ( : , 1:nJ )' * u;   e = RR( 1:nJ , 1:nJ ) \ uQ;
This computation of e can be done with nJ * ( 2 * m + nJ ) flops, if one takes
advantage of the triangular structure of RR( 1:nJ , 1:nJ ).

I do want to take advantage of the triangular stucture of RR. Since a typical
update of the QR factorization takes at most 6 * nJ * ( m + nJ ) flops, the QR
method should be much faster than the direct, e = A( : , J ) \ u, computation,
especially for large nJ.

My question is this. Should there be a speed advantage if I do
RR( 1:nJ , 1:nJ ) = matrix_type( RR( 1:nJ , 1:nJ ) , "upper" );
uQ = QQ( : , 1:nJ )' * u;   e = RR( 1:nJ , 1:nJ ) \ uQ;
on each iteration rather than doing just the last line? Or could I accomplish
the same thing by declaring
RR = matrix_type( RR , "upper" );
before the iterations start? Would this line force all subsequent  RR( 1:nJ ,
1:nJ )'s to be treated as upper triangular when e is computed?