Re: matrix_type problems

[David Bateman]

Vic Norton wrote:
> 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?
> 
> Thanks for your advice.
> 
> Regards,
> 
> Vic

Vic,

You don't have to force the recognition. The matrix will be recognized
as upper triangular. The slight performance advantage is only in the
fact that the code used to probe the matrix type in the \ operator
doesn't need to be called. So I'd expect a very minor speed up with the
use of the matrix_type function as the \ operator will treat RR as upper
triangular in any case. The same is not the case with versions of octave
earlier than 2.9.6.

Note also that the patch I sent means that "RR(1:nj, 1:nj) = ..."
invalidates the matrix type, so essentially RR(1:nj, 1:nj) =
matrix_type(RR(1:nj, 1:nj),'upper') is a null operation. Though RR =
matrix_type(RR,'upper') would make sense..

D.