lapackers

[James Demmel <demmel@cs.berkeley.edu>]

We didn't have a technical conversation, because the appropriate
technical folks at Mathworks were not on the call, but we agreed
to have such a call next week; we should try to schedule it so that
all interested parties can participate.

Also, Guillaume and Velvel spoke about the data, and found some
bugs that Guillaume is fixing and will rerun to collect new data.

We did have the following interesting "philosophical" discussion with
the Mathworks: Do we need to have just one QR decomposition,
that is also used in, say, the SVD and eigenroutines, where we don't care
about a nonnegative diagonal, if the alternative is (slightly) more 
accurate?
Or could we just have a new QR routine, with a new name,
that guarantees a nonnegative diagonal, rather than introducing
this change invisibly? Even though the new routine satisfies the
same specs as the old QR (which made no promises about the
diagonal), Mathworks' experience with their users is that any changes will
be noticed, and complained about, by someone. Of course such
rigorous adherence to backward compatibility can sometimes discourage
progress, but we could consider it (in this case!) if we can't
make slarfp at least as accurate as before.

Jim


Julien Langou wrote:
> Can we have a report on the phone call?
> (No problem if the answer is no.)
> --j
>
> On Fri, 7 May 2010, James Demmel wrote:
>
> > Ah, your point is that we get the cost benefit of scanning for zero
> > entries even if we use slarfg, and that we only introduced
> > scanning initially so that QR on an upper triangular input matrix
> > with some negative entries ran in O(n^2) time as expected,
> > not O(n^3).
> >
> > Still, the first stage of SVD on a tall skinny matrix is to call QR.
> > Then we'd need two flavors of QR (one calling slarfg and
> > one calling slarfp) in order for the QR in SVD to possibly
> > be slightly more accurate, and for other calls to QR to
> > return positive diagonals. We could do this.
> >
> > Jim
> >
> >
> >
> > Julien Langou wrote:
> > > Jim,
> > >
> > > I do not understand why you always associate
> > > 1) positive element on the diagonal of R, and,
> > > 2) scanning entries for narrow profile.
> > > This is two orthogonal contribution, isn't it?
> > >
> > > In any case, I trust your judgement for all this.
> > >
> > > ( Try to make a decision during the phone call !!! )
> > >
> > > Best,
> > > Julien.
> > >
> > >
> > >
> > > On Fri, 7 May 2010, James Demmel wrote:
> > >
> > > > Julien, thanks for the response.
> > > >
> > > > Your idea of keeping slarfp for QR alone, and dropping it in eig,
> > > > SVD, etc.,
> > > > is interesting. But I wonder about the case of the SVD of a tall 
> > > > skinny
> > > > matrix,
> > > > which starts by doing QR to reduce to a smaller matrix for bidiagonal
> > > > reduction.
> > > > It seems that we'd want the potential speed benefit of the new QR, if
> > > > the original
> > > > tall-skinny matrix also has a narrow profile (this could make the
> > > > overall SVD
> > > > much cheaper). I also wonder if we can easily separate which calls
> > > > should
> > > > be to the old slarfg and which to slarfp, it depends on the call trees
> > > > of all the
> > > > routines.
> > > >
> > > > So there is a set of design choices. We will be having a conference 
> > > > call
> > > > with the
> > > > Mathworks in about 30 minutes, and maybe get some feedback.
> > > >
> > > > Jim
> > > >
> > > >
> > > > Julien Langou wrote:
> > > > > Hi Jim,
> > > > >
> > > > > Thanks for the very clear update. I do not have much strong feelings.
> > > > >
> > > > > - I am in favor of "as simple as" possible codes. I did not consider
> > > > > the patch 1 necessarilly simple actually! But if patch 3 and 4 are
> > > > > even more complex ... having an extra routine to compute the norm
> > > > > diferently might be cleaner at the end.
> > > > >
> > > > > - I am not worried about a factor of 4x, 5x in loss of orthogonality
> > > > > additional to the O(n) anyhow.
> > > > >
> > > > > - We could also
> > > > >   - restore xlarfg in all the routines we do not know why we are 
> > > > > using
> > > > > xlarfp ... (eg: eig, svd, etc.) then for QR factorization I do not
> > > > > know maybe having the Alan Edleman QR factorization with xlarfp.
> > > > >
> > > > > best wishes, --julien.
> > > > >
> > > > >
> > > > >
> > > > >
> > > > > On Thu, 6 May 2010, James Demmel wrote:
> > > > >
> > > > > > I attach some interesting data Guillaume has collected from our work
> > > > > > on debugging slarfp. We are trying to come to a conclusion before
> > > > > > reporting progress to the Mathworks tomorrow. But depending on
> > > > > > our ultimate goal (cost vs accuracy), the data comparing various
> > > > > > possible
> > > > > > patches for slarfp is not conclusive, so we'd like your opinion.
> > > > > >
> > > > > > Briefly, there are several patches of different costs and worst-case
> > > > > > accuracies that are much better than the buggy version, and probably
> > > > > > "good enough". But all are less accurate than the original 
> > > > > > slarfg, by
> > > > > > about 4x or 5x (eyeball metric from the plots; more exact data is
> > > > > > available).
> > > > > > Whether this lower accuracy is worth worrying about is a judgment
> > > > > > call about which I'd like your opinion.
> > > > > >
> > > > > > To remind you of notation, the input is the vector [alpha;x] and the
> > > > > > output is the orthogonal matrix Q = I - tau*[1;v]*[1,v'] and beta 
> > > > > > >=0
> > > > > > such that Q*[alpha;x] = [beta;0].
> > > > > >
> > > > > > The attached pdf has 3 pages, showing histograms of errors over
> > > > > > certain "torture test" suites. The error shown is loss of
> > > > > > orthogonality,
> > > > > > measured as norm(Q'*Q-I), computed in double. (There is also
> > > > > > data for the residual norm(Q*[alpha;x] - [beta;0])/beta, but
> > > > > > this is enough data for now). More precisely, we histogram
> > > > > >   log10(error/macheps)
> > > > > > i.e. how many decimal digits we lose over the "best" error of
> > > > > > macheps = 2^(-24). Actually, a large number of test cases had
> > > > > > smaller errors than this, and are not shown, but counted as
> > > > > > "not shown" in the caption to each figure.
> > > > > >
> > > > > > There are 3 pages of plots corresponding to 3 kinds of tests.
> > > > > > The input vector is of the form [alpha;x], and the dimension
> > > > > > of x can be 2, 100 and 1000. When the dimension is 2,
> > > > > > each component alpha and x(i) is chosen from a large set
> > > > > > of floating point numbers spanning the range from 0,
> > > > > > through the denorms, and from underflow to overflow threshholds,
> > > > > > with various mantissas and signs. When the dimension
> > > > > > of x is 100 or 1000, alpha and x(1) are chosen similarly generally,
> > > > > > but x(2:end) are chosen very tiny compared to x(1) to maximize
> > > > > > the relative error in computing snrm2(x) (i.e. x(2:end) are all
> > > > > > about sqrt(macheps)*x(1)).
> > > > > >
> > > > > > For each page of tests, 5 plots are shown:
> > > > > >   (a) slarfg: original implementation, before our "improvements"
> > > > > >   (b) slarfp (3.2, buggy): the version from LAPACK 3.2 that
> > > > > >         caused the bug report
> > > > > >   (c) slarfp (patch 1): the version that Julien and I worked on
> > > > > >   (d) slarfp (patch 2): not a practical patch, rather it uses a 
> > > > > > bit of
> > > > > >             double precision in a subset of the code that we
> > > > > > identified
> > > > > >             as problematic from our delta-debugging work
> > > > > >   (e) slarfp (patch 3): a suggestion of Velvel, where the output tau
> > > > > >             is explicitly computed from its definition
> > > > > > 2/(1+norm(v,2)^2)
> > > > > >             and so costing an additional O(n) work
> > > > > >   (f) slarfp (patch 4): a suggestion of Guillaume, that also
> > > > > > recomputes
> > > > > >             a vector norm, but differently.
> > > > > >
> > > > > > The main features of these plots are as follows:
> > > > > >    The bug in slarfp  (3.2) shows up clearly.
> > > > > >    For dimension n = 100 and n=1000, the worst case error in slarfg
> > > > > >       and all patched versions of slarfp grows like O(n*macheps) as
> > > > > > expected,
> > > > > >       since this is the error from snrm2. We have lived with this
> > > > > > situation
> > > > > >       (in slarfg) for a long time, and it dominates the differences
> > > > > > between the
> > > > > >       different patches. We could decide that because of this, the
> > > > > > difference
> > > > > >       between the patches is too small to matter, and we should
> > > > > > pick the
> > > > > >       simplest and cheapest (patch 1).
> > > > > >    For dimension n=2, the worst case error for slarfg and slarfp
> > > > > > (patch 3)
> > > > > >       are about the same, and about 4x smaller than slarfp 
> > > > > > (patches 1
> > > > > > and 2)
> > > > > >       and 20x smaller than patch 4. (The 4x and 20x are "eyeball"
> > > > > > values, Guillaume
> > > > > >       has the exact data.)
> > > > > >   For dimension n=100, slarfg is about 5x (eyeball) better than
> > > > > >       patches 1, 2 and 3, which are about the same, but only 2x
> > > > > > better than
> > > > > >       patch 4.
> > > > > >   For dimension n = 1000, the situation is similar, but now patch 
> > > > > > 4 is
> > > > > > about as good
> > > > > >      as slarfg.
> > > > > >
> > > > > > So to repeat my question: do we care that the simplest patch, and
> > > > > > sometimes the
> > > > > > fancier ones, are as much as 4x to 5x worse than the original slarfg
> > > > > > (epending
> > > > > > on dimension)?
> > > > > >
> > > > > > Jim