[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: QR updating changeset
From: |
Jaroslav Hajek |
Subject: |
Re: QR updating changeset |
Date: |
Wed, 5 Mar 2008 11:24:25 +0100 |
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