To avoid having the strategy method being overriden in the constructor, we
could instead implement a constructor that takes a column Comparator that
determines if two columns should be exchanged at each stage.
In the interest of maintaining a clean interface, I don't know if this
constructor should be public though.
Chris
On 7 August 2011 16:43, Greg Sterijevski <gsterijev...@gmail.com> wrote:
> Chris,
>
> In regard to the pivoting, do you think that it would be useful to allow
> subclasses to pivot using other strategies? The pivoting I have looks for
> the next column with the largest norm. For most garden variety problems
> this
> will be okay. However, you can (I am not sure how effective this will be)
> choose the column which has the lowest [average] inner product (ie is least
> correlated in some sense to the remaining columns).
>
> The easiest way to accomplish this is to present the remainder columns to
> some method, and have that method return an index... The method would be
> protected so that it could be overridden by classes wishing to modify that
> behavior. The only problem I see is that the actual decomposition would
> need
> to be moved out of the constructor (we would have an overridable method
> being called in the constructor).
>
> -Greg
>
> On Sun, Aug 7, 2011 at 7:29 AM, Chris Nix <chris....@gmail.com> wrote:
>
> > Oooops, the ********* below should read MATH-630.
> >
> > Chris
> >
> > On 7 August 2011 13:28, Chris Nix <chris....@gmail.com> wrote:
> >
> > > Thanks, Greg, for looking more at this. Apologies I've not been able
> to
> > > focus on this too much recently.
> > >
> > > Yes, the approach above works, however the solvers also require a
> change
> > so
> > > that they don't unexpectedly solve for a permuted matrix.
> Additionally,
> > a
> > > getRank() method can now be implemented much more efficiently than the
> > > getRank() from SingularValueDecomposition.
> > >
> > > I submitted an initial patch at *********** with these bits working,
> > > however this patch introduces a new RRQRDecomposition interface
> extending
> > > QRDecomposition. In light of the insights above from the team, I'll
> > > implement it instead within the existing class structure and re-submit.
> > >
> > > If the pivoting code is to be included in QRDecomposition, then perhaps
> > an
> > > extra constructor is required to perform column pivoting since the RRQR
> > > decomposition of matrix A produces Q, R and P such A = QRP^T and would
> > break
> > > any existing code that expects a plain decomposition of A=QR.
> > >
> > > Chris.
> > >
> > > On 6 August 2011 21:34, Greg Sterijevski <gsterijev...@gmail.com>
> wrote:
> > >
> > >> Hello All,
> > >>
> > >> Not sure if this is stepping on toes, but here is what I thought of to
> > >> deal
> > >> with pivoting. I would only need to modify the constructor of the
> > >> QRDecompImpl class:
> > >>
> > >> public QRDecompositionPivotImpl(RealMatrix matrix) {
> > >>
> > >> final int m = matrix.getRowDimension();
> > >> final int n = matrix.getColumnDimension();
> > >> final int mn = FastMath.min(m, n);
> > >>
> > >> qrt = matrix.transpose().getData();
> > >> rDiag = new double[FastMath.min(m, n)];
> > >> cachedQ = null;
> > >> cachedQT = null;
> > >> cachedR = null;
> > >> cachedH = null;
> > >>
> > >>
> > >> /*
> > >> * The QR decomposition of a matrix A is calculated using
> > >> Householder
> > >> * reflectors by repeating the following operations to each
> minor
> > >> * A(minor,minor) of A:
> > >> */
> > >>
> > >> pivots = new int[n];
> > >> for (int i = 0; i < n; i++) {
> > >> pivots[i] = i;
> > >> }
> > >>
> > >> int pivotIdx = -1;
> > >> double bestNorm = 0.0;
> > >> for (int minor = 0; minor < mn; minor++) {
> > >> bestNorm = 0.0;
> > >> pivotIdx = 0;
> > >> for (int i = minor; i < n; i++) {
> > >> final double[] qrtMinor = qrt[i];
> > >> double xNormSqr = 0.0;
> > >> for (int row = minor; row < m; row++) {
> > >> final double c = qrtMinor[row];
> > >> xNormSqr += c * c;
> > >> }
> > >> if (xNormSqr > bestNorm) {
> > >> bestNorm = xNormSqr;
> > >> pivotIdx = i;
> > >> }
> > >> }
> > >>
> > >>
> > >> if( pivotIdx != minor){
> > >> int pvt = pivots[minor];
> > >> pivots[minor] = pivots[pivotIdx];
> > >> double[] qswp = qrt[minor];
> > >> qrt[minor] = qrt[pivotIdx];
> > >> for( int i = minor + 1; i < n; i++){
> > >> if( i <= pivotIdx ){
> > >> final int itmp = pivots[i];
> > >> pivots[i] = pvt;
> > >> pvt = itmp;
> > >> final double[] tmp = qrt[i];
> > >> qrt[i] = qswp;
> > >> qswp=tmp;
> > >> }
> > >> }
> > >> }
> > >>
> > >> final double[] qrtMinor = qrt[minor];
> > >> /*
> > >> * Let x be the first column of the minor, and a^2 = |x|^2.
> > >> * x will be in the positions qr[minor][minor] through
> > >> qr[m][minor].
> > >> * The first column of the transformed minor will be
> > >> (a,0,0,..)'
> > >> * The sign of a is chosen to be opposite to the sign of
> the
> > >> first
> > >> * component of x. Let's find a:
> > >> */
> > >> final double a = (qrtMinor[minor] > 0) ?
> > >> -FastMath.sqrt(bestNorm) : FastMath.sqrt(bestNorm);
> > >> rDiag[minor] = a;
> > >>
> > >> if (a != 0.0) {
> > >>
> > >> /*
> > >> * Calculate the normalized reflection vector v and
> > >> transform
> > >> * the first column. We know the norm of v beforehand:
> v
> > =
> > >> x-ae
> > >> * so |v|^2 = <x-ae,x-ae> = <x,x>-2a<x,e>+a^2<e,e> =
> > >> * a^2+a^2-2a<x,e> = 2a*(a - <x,e>).
> > >> * Here <x, e> is now qr[minor][minor].
> > >> * v = x-ae is stored in the column at qr:
> > >> */
> > >> qrtMinor[minor] -= a; // now |v|^2 =
> > -2a*(qr[minor][minor])
> > >>
> > >> /*
> > >> * Transform the rest of the columns of the minor:
> > >> * They will be transformed by the matrix H =
> > I-2vv'/|v|^2.
> > >> * If x is a column vector of the minor, then
> > >> * Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x -
> > 2<x,v>/|v|^2
> > >> v.
> > >> * Therefore the transformation is easily calculated by
> > >> * subtracting the column vector (2<x,v>/|v|^2)v from
> x.
> > >> *
> > >> * Let 2<x,v>/|v|^2 = alpha. From above we have
> > >> * |v|^2 = -2a*(qr[minor][minor]), so
> > >> * alpha = -<x,v>/(a*qr[minor][minor])
> > >> */
> > >> for (int col = minor + 1; col < n; col++) {
> > >> final double[] qrtCol = qrt[col];
> > >> double alpha = 0;
> > >> for (int row = minor; row < m; row++) {
> > >> alpha -= qrtCol[row] * qrtMinor[row];
> > >> }
> > >> alpha /= a * qrtMinor[minor];
> > >>
> > >> // Subtract the column vector alpha*v from x.
> > >> for (int row = minor; row < m; row++) {
> > >> qrtCol[row] -= alpha * qrtMinor[row];
> > >> }
> > >> }
> > >> }
> > >> }
> > >> }
> > >>
> > >>
> > >> All I am doing is looking forward and taking the next column with the
> > >> largest norm. Then I rearrange the Q's. Is this what you had in mind
> > >> Chris?
> > >> The result will be Q,R and the pivot array for which we can implement
> a
> > >> getter?
> > >>
> > >> -Greg
> > >>
> > >
> > >
> >
>
