Re: [eigen] Golub-Kahan bidiagonalization

["Schleimer, Ben" <bensch128@xxxxxxxxx>]

Hi Benoit,
   It might be useful for there to be a third bidiagonalization class which takes in an epsilon. It would run the partial pivoting version first and then check the error. If it's greater then epsilon, then it would run the full pivoting solution. I don't know how much overhead the error calculation would add but it might be worth it to make the api simpler..

just my two cents,
Ben



--- On Tue, 4/28/09, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:

> From: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>
> Subject: [eigen] Golub-Kahan bidiagonalization
> To: eigen@xxxxxxxxxxxxxxxxxxx
> Date: Tuesday, April 28, 2009, 6:28 AM
> Hi,
> 
> this is mostly a followup to an old conversation with Keir
> as things
> are a little clearer in my mind now, but also to gather
> some feedback
> if there is an expert around.
> 
> what i discuss here applies equally to tridiagonalization
> and
> Hessenberg reduction.
> 
> 1) instability in the Golub-Kahan bidiagonalization
> algorithm.
> 
> This is algorithm 5.4.2 in Golub&vanLoan, called
> "Householder
> bidiagonalization".
> 
> It consists in left-multiplying by a Householder reflection
> to put the
> first column in the form (x,0,...,0), then
> right-multiplying by
> another Householder reflection to put the remainder of the
> first row
> (that is, without the 1st entry) in the form (y,0,....,0),
> and so on.
> 
> The problem that I see is that if some column is close to
> zero, then a
> small perturbation of it results in a large perturbation of
> the
> corresponding Householder vector. This is because the
> Householder
> vector is normalized, and so this amounts to normalizing a
> near-zero
> vector. Additional instability comes from the fact that
> this involves
> taking the square root of a near-zero number (to compute
> the norm).
> 
> 2) proposed solution: full pivoting
> 
> Alongside plain bidiagonalization, we could offer a variant
> doing as
> follows. It would be slower but more reliable.
> 
> At every step, before computing the Householder vector of a
> column, we
> could look up the (squared) norms of all remaining columns
> and permute
> columns to make the biggest remaining column the one that
> we are
> putting in form (x,0....,0). Same for rows. These
> permutation matrices
> being unitaries, fit perfectly in the decomposition which
> is already
> of the form UAV with U,V two independent unitaries.
> 
> Since Golub-Kahan bidiagonalization already has (8/3)n^3
> complexity,
> the overhead of full pivoting here would be a bit lower
> than it is for
> LU.
> 
> 3) API for that.
> 
> standard Bidiagonalization and pivoting Bidiagonalization
> would be 2
> separate classes, but they could have a common base class
> abstracting
> that difference away from the SVD. The QR iteration to
> compute the SVD
> is the same in both cases.
> 
> Cheers,
> Benoit
> 
> 
>