QR updating changeset

[John W. Eaton]

On 24-Feb-2008, Jaroslav Hajek wrote:

| hello,
| 
| attached is a changeset to implement updating QR factorizations in Octave.
| Brief outline: adds new Fortran sources library qrupdate in libcruft/,
| implements updating methods for the QR and ComplexQR classes,
| and adds DEFUNs qrupdate, qrinsert and qrdelete.
| Each of the DEFUNs contains basic tests, but more detailed testing of corner
| cases is certainly desirable (e.g., I'm not quite sure how the functions 
should
| behave give zero-sized args).
| 
| The behaviour is slightly more extended over the Matlab's, namely that 
updating
| and inserting columns do not gripe about non-square basis matrix, but proceeds
| with the projection of the given vector onto the basis instead (this
| comes out as
| a natural extension in the algorithm, and is IMHO an useful extension).
| 
| another quirk is qrdelete with 'col' vs. 'col+' - see documentation.

I applied this patch with some changes but would like for you to make
some more based on my comments below.

First, I put qrupdate, qrinsert, qrdelete in the qr.cc file instead of
in their own files.

I changed int to octave_idx_type in several places.

I changed

| +              retval (0) = fact.Q ();
| +              retval (1) = fact.R ();

to

              retval(1) = fact.R ();
              retval(0) = fact.Q ();

so that retval only has to be resized once.

I rewrote

| +void QR::economize ()
| +{
| +  idx_vector c (':'), i (Range (1, r.columns ()));
| +  q = Matrix (q.index (c, i));
| +  r = Matrix (r.index (i, c));
| +}

(and the ComplexQR version) like this:

  void
  QR::economize (void)
  {
    octave_idx_type r_nc = r.columns ();

    if (r.rows () > r_nc)
      {
        q.resize (q.rows (), r_nc);
        r.resize (r_nc, r_nc);
      }
  }

The resize method should be more efficient than using the index
function.

I added the check on the size of R since it seems that resizing would
cause trouble if R has fewer rows than columns, since then I think you
would be resizing Q to have more columns that it originally has (in
your version, you'd be indexing outside the bounds of Q).  Or am I
missing something?

In the qrupdate, qrinsert, and qrdelete functions, you have code like
this to deal with the arguments:

| +  octave_value argQ,argR,argj,argor;
| +
| +  if ((nargin == 3 || nargin == 4) && nargout == 2
| +      && (argQ = args (0), argQ.is_matrix_type ())
| +      && (argR = args (1), argR.is_matrix_type ())
| +      && (argj = args (2), argj.is_scalar_type ())
| +      && (nargin < 4 || (argor = args (3), argor.is_string ())))
| +    {
| +      int m = argQ.rows (), k = argQ.columns (), n = argR.columns ();
| +      bool row = false, colp = false;
| +      std::string orient = (nargin < 4) ? "col" : argor.string_value ();
| +
| +      if (nargin < 4 || (row = orient == "row") 
| +          || (orient == "col") || (colp = orient == "col+"))
| +        if (argR.rows() == k)

Will you please fix the above to avoid using comma-expressions and
performing assignments inside the if conditions?  Also, we typically
try to avoid mixed-case variables names in the Octave sources.  We
also don't usually check nargout unless it is to provide different
output values depending on how many were requested.

Finally, I'm seeing this error now (even without my changes, so I
don't think it is a problem that I introduced):

   ====>>>>> processing /tmp/jwe/octave/src/DLD-FUNCTIONS/qrinsert.cc
    ***** test
   A = [0.620405 + 0.956953i  0.480013 + 0.048806i  0.402627 + 0.338171i;
        0.589077 + 0.658457i  0.013205 + 0.279323i  0.229284 + 0.721929i;
        0.092758 + 0.345687i  0.928679 + 0.241052i  0.764536 + 0.832406i;
        0.912098 + 0.721024i  0.049018 + 0.269452i  0.730029 + 0.796517i;
        0.112849 + 0.603871i  0.486352 + 0.142337i  0.355646 + 0.151496i ];

   x = [0.20351 + 0.05401i;
        0.13141 + 0.43708i;
        0.29808 + 0.08789i;
        0.69821 + 0.38844i;
        0.74871 + 0.25821i ];

   [Q,R] = qr(A);
   [Q,R] = qrinsert(Q,R,3,x);
   assert(norm(vec(Q'*Q - eye(5)),Inf) < 1e1*eps)
   assert(norm(vec(triu(R)-R),Inf) == 0)
   assert(norm(vec(Q*R - [A(:,1:2) x A(:,3)]),Inf) < norm(A)*1e1*eps)

  !!!!! test failed
  error: assert (norm (vec (Q * R - [A(:, 1:2), x, A(:, 3)]), Inf) < norm (A) * 
1e1 * eps) failed

Your patch and my changes are all in my public archive now, so will
you please update and take a look?

Thanks,

jwe