Re: QR updating changeset

[Jaroslav Hajek]

On Wed, Mar 5, 2008 at 3:52 AM, John W. Eaton <address@hidden> wrote:
> 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.

Should this also be done when the integer argument doesn't have a
meaning of an index or size? (e.g., a LAPACK info argument)


>
> 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.
>

Nice, I'll remember that.

> 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.
>

Thanks, I didn't realize that resize can preserve data.

> 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?
>

No, the correction is right.

> 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.
>

OK, I'll revamp this.
John, I know it would be a boring work, but could you please consider
summarizing coding style recommendations for Octave somewhere? (e.g.
README.devel). For instance, start with a reference to GNU coding
standards and then summarize whatever goes beyond those (I don't
recall anything about comma exps there, so that's a good example).

Alternatively, I can start this job by gathering all the guidelines
I've already learned from you, and leave it to you to add and correct
as necessary afterwards.

>
> 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?
>

Strange, I think I recall all tests pass. I'll try to reproduce and
fix the error.

> Thanks,
>
> jwe
>
>
>

regards,

-- 
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz