You are viewing a plain text version of this content. The canonical link for it is here.
Posted to dev@commons.apache.org by Greg Sterijevski <gs...@gmail.com> on 2011/08/06 22:34:12 UTC
Re: [math] Pivoting in QR decomposition
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
Re: [math] Pivoting in QR decomposition
Posted by Greg Sterijevski <gs...@gmail.com>.
Your selection methodology might be more complicated than pairwise
comparison. So I think you would have to pass the entire set of remaining
columns. That is an implementation detail, however.
On the interface, yes the public face should not expose the constructor
which takes the PivotStrategyComparator. It should be part of a protected
constructor or, I suppose, we could have a static factory class which spits
out different implementations.
-Greg
On Sun, Aug 7, 2011 at 11:33 AM, Chris Nix <ch...@gmail.com> wrote:
> 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 <gs...@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 <ch...@gmail.com> wrote:
> >
> > > Oooops, the ********* below should read MATH-630.
> > >
> > > Chris
> > >
> > > On 7 August 2011 13:28, Chris Nix <ch...@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 <gs...@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
> > > >>
> > > >
> > > >
> > >
> >
>
Re: [math] Pivoting in QR decomposition
Posted by Chris Nix <ch...@gmail.com>.
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 <gs...@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 <ch...@gmail.com> wrote:
>
> > Oooops, the ********* below should read MATH-630.
> >
> > Chris
> >
> > On 7 August 2011 13:28, Chris Nix <ch...@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 <gs...@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
> > >>
> > >
> > >
> >
>
Re: [math] Pivoting in QR decomposition
Posted by Greg Sterijevski <gs...@gmail.com>.
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 <ch...@gmail.com> wrote:
> Oooops, the ********* below should read MATH-630.
>
> Chris
>
> On 7 August 2011 13:28, Chris Nix <ch...@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 <gs...@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
> >>
> >
> >
>
Re: [math] Pivoting in QR decomposition
Posted by Chris Nix <ch...@gmail.com>.
Oooops, the ********* below should read MATH-630.
Chris
On 7 August 2011 13:28, Chris Nix <ch...@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 <gs...@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
>>
>
>
Re: [math] Pivoting in QR decomposition
Posted by Chris Nix <ch...@gmail.com>.
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 <gs...@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
>