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/07/21 06:27:12 UTC
[math] Pivoting in QR decomposition
Hello,
In attempting to build constrained ols into the OLSMultipleLinearRegression
class, I rediscovered the fact that to keep the regression coefficients in
their canonical order, I need to know how the QR decomposition pivoted. In
other words, while I get the correct constrained parameters, their order is
changed. Is there any way to get this info out (pivot order)?
Thanks,
-Greg
Re: [math] Pivoting in QR decomposition
Posted by Ted Dunning <te...@gmail.com>.
Sounds like you need a new method.
It also sounds like the Q and R components of the decomposition (if
accessible) will be rendered in permuted order. Is there such an accessor?
If so, is the doc up to snuff on this caveat?
On Wed, Jul 20, 2011 at 9:27 PM, Greg Sterijevski <gs...@gmail.com>wrote:
> Is there any way to get this info out (pivot order)?
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
>
Re: [math] Pivoting in QR decomposition
Posted by Greg Sterijevski <gs...@gmail.com>.
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>.
I second +1.
On Sat, Jul 23, 2011 at 4:06 PM, Phil Steitz <ph...@gmail.com> wrote:
> On 7/23/11 1:37 PM, Ted Dunning wrote:
> > Also, if the QRDecomposition interface *is* extended with getP, it is the
> > work of a moment to change it to an abstract class instead of an
> interface
> > and provide a default method for getP that throws
> > UnsupportedOperationException. Since there are very unlikely to be any
> user
> > extensions of QRDecomposition in the wild, we might as well do both
> changes
> > at once. For that matter, getP can also easily return an identity
> > permutation by default.
> >
> > Much like an immutable set or list implements the standard remove method,
> > but complains when it is called, it is quite reasonable for
> QRDecomposition
> > to not be quite the lowest common denominator, but something more like a
> > mid-line that contains the things that you might expect a QR
> decomposition
> > to do. Having a few (very few) methods of this sort is preferable to
> having
> > an elaborate class structure of alternative interfaces.
>
> +1
>
> Looks like a reasonable approach to me, with getP returning identity
> by default.
>
> Phil
> >
> > On Sat, Jul 23, 2011 at 12:45 PM, Chris Nix <ch...@gmail.com> wrote:
> >
> >> We can do this with the current Householder reflection implementation,
> >> except instead of just obtaining reflections from columns in sequence
> >> across
> >> the input matrix, we select the column with the greatest L2-norm at each
> >> iteration. The resulting permutation matrix is thus built and can be
> >> returned later with a getP() method. A pleasing by-product is that the
> >> resulting R matrix is 'rank-revealing', allowing for a quicker getRank()
> >> than currently exists in the SingularValueDecompositionImpl class.
> >>
> >> Does it make sense to extend the current QRDecomposition interface to
> one
> >> for rank-revealing QR decompositions that have the existing methods,
> plus a
> >> getP() and getRank()? The implementing class could extend the current
> >> QRDecompositionImpl class to save reproducing code, at the cost of
> opening
> >> up some private variables and methods, and the solver.
> >>
> >> I'll open a JIRA issue.
> >>
> >> Chris.
> >>
> >> On 23 July 2011 18:44, Greg Sterijevski <gs...@gmail.com> wrote:
> >>
> >>> Chris, you had an algorithm in mind? -Greg
> >>>
> >>> On Sat, Jul 23, 2011 at 11:29 AM, Phil Steitz <ph...@gmail.com>
> >>> wrote:
> >>>
> >>>> On 7/22/11 11:40 AM, Greg Sterijevski wrote:
> >>>>> Sorry,
> >>>>>
> >>>>> public ConstrainedOLSMultipleRegression extends
> >> OLSMultipleRegression{}
> >>>>> should read:
> >>>>>
> >>>>> public ConstrainedOLSMultipleRegression extends
> >> OLSMultipleRegression{
> >>>>> @Override
> >>>>> public void newSampleData(double[] data, double[][] coeff,
> >> double[]
> >>>> rhs,
> >>>>> int nob, int nvars) {
> >>>>> adjustData( data, coeff, rhs);
> >>>>> super.newSampleData(data, nobs, nvars);
> >>>>> qr = new QRDecompositionImpl(X);
> >>>>> }
> >>>>>
> >>>>>> }
> >>>>> The data would be transformed on the way in, and everything else
> >> would
> >>>>> remain the same...
> >>>> Got it. Sounds good. Patch away...
> >>>>
> >>>> Couple of things to keep in mind:
> >>>>
> >>>> 0) We may want to dispense with the QRDecomposition interface
> >>>> altogther. If we keep it, we should trim it down to what is common
> >>>> and meaningfully implemented in all impls. So both the Householder
> >>>> and permutation getters are dropped. If you need a pivoting impl,
> >>>> go a head and code it up and we can reassess the interface.
> >>>>
> >>>> 1) We should be aiming to standardize the regression API. Lets pick
> >>>> up the other thread on regression refactoring. Before hacking too
> >>>> much more on OLS, lets refine and retrofit the new regression API on
> >>>> this class.
> >>>>
> >>>> Phil
> >>>>>
> >>>>> On Fri, Jul 22, 2011 at 1:37 PM, Greg Sterijevski <
> >>>> gsterijevski@gmail.com>wrote:
> >>>>>> On the need for pivoting:
> >>>>>>
> >>>>>> Here is my first approach for changing OLSMultipleRegression to do
> >>>>>> constrained estimation:
> >>>>>>
> >>>>>> public double[] calculateBeta(double[][] coeff, double[] rhs) {
> >>>>>> if (rhs.length != coeff.length) {
> >>>>>> throw new IllegalArgumentException("");
> >>>>>> }
> >>>>>> for (double[] rest : coeff) {
> >>>>>> if (rest.length != this.X.getColumnDimension()) {
> >>>>>> throw new IllegalArgumentException("");
> >>>>>> }
> >>>>>> }
> >>>>>> RealMatrix Coeff = new Array2DRowRealMatrix(coeff, false);
> >>>>>> RealVector rhsVec = new ArrayRealVector(rhs);
> >>>>>> QRDecomposition coeffQRd = new
> >>>>>> QRDecompositionImpl(Coeff.transpose());
> >>>>>> RealMatrix Qcoeff = coeffQRd.getQ();
> >>>>>> RealMatrix R = X.multiply(Qcoeff);
> >>>>>>
> >>>>>> final int nvars = X.getColumnDimension();
> >>>>>> final int nobs = X.getRowDimension();
> >>>>>> final int ncons = coeff.length;
> >>>>>>
> >>>>>> RealMatrix R2 = R.getSubMatrix(
> >>>>>> 0, nobs - 1, ncons, nvars - 1);
> >>>>>>
> >>>>>> RealMatrix R1 = R.getSubMatrix(
> >>>>>> 0, nobs - 1, 0, ncons - 1);
> >>>>>>
> >>>>>> RealVector gamma = rhsVec.copy();
> >>>>>>
> >>>>>> RealMatrix coeffR = coeffQRd.getR().getSubMatrix(
> >>>>>> 0, ncons - 1, 0, ncons - 1);
> >>>>>>
> >>>>>> MatrixUtils.solveLowerTriangularSystem(coeffR.transpose(),
> >>>> gamma);
> >>>>>> RealVector gammPrime = Y.subtract(R1.operate(gamma));
> >>>>>>
> >>>>>> QRDecomposition qr2 = new QRDecompositionImpl(R2);
> >>>>>>
> >>>>>> RealVector constrainedSolution =
> >>>>>> (qr2.getSolver().solve(gammPrime));
> >>>>>>
> >>>>>> RealVector stackedVector =
> >>>>>> new ArrayRealVector(
> >>>>>> gamma.toArray(),
> >>>>>> constrainedSolution.toArray());
> >>>>>>
> >>>>>> stackedVector = Qcoeff.operate(stackedVector);
> >>>>>>
> >>>>>> return stackedVector.toArray();
> >>>>>> }
> >>>>>>
> >>>>>> This approach is based on Dongarra et al:
> >>>>>>
> >>>>>> LAPACK Working Note
> >>>>>> Generalized QR Factorization and its Applications
> >>>>>> Work in Progress
> >>>>>> E. Anderson, Z. Bai and J. Dongarra
> >>>>>> December 9, 1991
> >>>>>> August 9, 1994
> >>>>>>
> >>>>>> There is nothing terrible about this approach, the coding is not
> >>>> finished
> >>>>>> and tidy, but its a work in progress.
> >>>>>>
> >>>>>> I am also aware of second approach. I do not have a cite for it, I
> >>> think
> >>>> I
> >>>>>> may have derived it myself, but it would not surprise me if it is in
> >>>> some
> >>>>>> textbook somewhere... That second approach takes the QR
> >> decomposition
> >>> of
> >>>> the
> >>>>>> coefficient matrix and calculates adjustment matrices for the design
> >>>> matrix
> >>>>>> and dependent vector. The problem is that I need to reorganize the
> >>>> design
> >>>>>> matrix by the pivots of the QR decomposition. Once I have the
> >>> adjustment
> >>>>>> matrices, everything should proceed as in the case of an
> >> unconstrained
> >>>>>> estimation. I like the idea that if we transform the data,
> >> everything
> >>>> works
> >>>>>> the same way.
> >>>>>>
> >>>>>> Since then the ConstrainedOLSMultipleRegression class looks like:
> >>>>>> public ConstrainedOLSMultipleRegression extends
> >> OLSMultipleRegression{
> >>>>>> }
> >>>>>>
> >>>>>>
> >>>>>> As for the fact that the QRDecompositionImpl reflects its interface.
> >>> We
> >>>>>> should probably add the functions:
> >>>>>> public int[] getPivots();
> >>>>>> public boolean isPivotting();
> >>>>>>
> >>>>>> to the interface. As Christopher pointed out, if the current
> >>>> decomposition
> >>>>>> is non pivoting, its pivot record is the canonical one,
> >>> {0,1,2,...,n-1}.
> >>>>>> -Greg
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>
> >>>> ---------------------------------------------------------------------
> >>>> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> >>>> For additional commands, e-mail: dev-help@commons.apache.org
> >>>>
> >>>>
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
>
Re: [math] Pivoting in QR decomposition
Posted by Phil Steitz <ph...@gmail.com>.
On 7/23/11 1:37 PM, Ted Dunning wrote:
> Also, if the QRDecomposition interface *is* extended with getP, it is the
> work of a moment to change it to an abstract class instead of an interface
> and provide a default method for getP that throws
> UnsupportedOperationException. Since there are very unlikely to be any user
> extensions of QRDecomposition in the wild, we might as well do both changes
> at once. For that matter, getP can also easily return an identity
> permutation by default.
>
> Much like an immutable set or list implements the standard remove method,
> but complains when it is called, it is quite reasonable for QRDecomposition
> to not be quite the lowest common denominator, but something more like a
> mid-line that contains the things that you might expect a QR decomposition
> to do. Having a few (very few) methods of this sort is preferable to having
> an elaborate class structure of alternative interfaces.
+1
Looks like a reasonable approach to me, with getP returning identity
by default.
Phil
>
> On Sat, Jul 23, 2011 at 12:45 PM, Chris Nix <ch...@gmail.com> wrote:
>
>> We can do this with the current Householder reflection implementation,
>> except instead of just obtaining reflections from columns in sequence
>> across
>> the input matrix, we select the column with the greatest L2-norm at each
>> iteration. The resulting permutation matrix is thus built and can be
>> returned later with a getP() method. A pleasing by-product is that the
>> resulting R matrix is 'rank-revealing', allowing for a quicker getRank()
>> than currently exists in the SingularValueDecompositionImpl class.
>>
>> Does it make sense to extend the current QRDecomposition interface to one
>> for rank-revealing QR decompositions that have the existing methods, plus a
>> getP() and getRank()? The implementing class could extend the current
>> QRDecompositionImpl class to save reproducing code, at the cost of opening
>> up some private variables and methods, and the solver.
>>
>> I'll open a JIRA issue.
>>
>> Chris.
>>
>> On 23 July 2011 18:44, Greg Sterijevski <gs...@gmail.com> wrote:
>>
>>> Chris, you had an algorithm in mind? -Greg
>>>
>>> On Sat, Jul 23, 2011 at 11:29 AM, Phil Steitz <ph...@gmail.com>
>>> wrote:
>>>
>>>> On 7/22/11 11:40 AM, Greg Sterijevski wrote:
>>>>> Sorry,
>>>>>
>>>>> public ConstrainedOLSMultipleRegression extends
>> OLSMultipleRegression{}
>>>>> should read:
>>>>>
>>>>> public ConstrainedOLSMultipleRegression extends
>> OLSMultipleRegression{
>>>>> @Override
>>>>> public void newSampleData(double[] data, double[][] coeff,
>> double[]
>>>> rhs,
>>>>> int nob, int nvars) {
>>>>> adjustData( data, coeff, rhs);
>>>>> super.newSampleData(data, nobs, nvars);
>>>>> qr = new QRDecompositionImpl(X);
>>>>> }
>>>>>
>>>>>> }
>>>>> The data would be transformed on the way in, and everything else
>> would
>>>>> remain the same...
>>>> Got it. Sounds good. Patch away...
>>>>
>>>> Couple of things to keep in mind:
>>>>
>>>> 0) We may want to dispense with the QRDecomposition interface
>>>> altogther. If we keep it, we should trim it down to what is common
>>>> and meaningfully implemented in all impls. So both the Householder
>>>> and permutation getters are dropped. If you need a pivoting impl,
>>>> go a head and code it up and we can reassess the interface.
>>>>
>>>> 1) We should be aiming to standardize the regression API. Lets pick
>>>> up the other thread on regression refactoring. Before hacking too
>>>> much more on OLS, lets refine and retrofit the new regression API on
>>>> this class.
>>>>
>>>> Phil
>>>>>
>>>>> On Fri, Jul 22, 2011 at 1:37 PM, Greg Sterijevski <
>>>> gsterijevski@gmail.com>wrote:
>>>>>> On the need for pivoting:
>>>>>>
>>>>>> Here is my first approach for changing OLSMultipleRegression to do
>>>>>> constrained estimation:
>>>>>>
>>>>>> public double[] calculateBeta(double[][] coeff, double[] rhs) {
>>>>>> if (rhs.length != coeff.length) {
>>>>>> throw new IllegalArgumentException("");
>>>>>> }
>>>>>> for (double[] rest : coeff) {
>>>>>> if (rest.length != this.X.getColumnDimension()) {
>>>>>> throw new IllegalArgumentException("");
>>>>>> }
>>>>>> }
>>>>>> RealMatrix Coeff = new Array2DRowRealMatrix(coeff, false);
>>>>>> RealVector rhsVec = new ArrayRealVector(rhs);
>>>>>> QRDecomposition coeffQRd = new
>>>>>> QRDecompositionImpl(Coeff.transpose());
>>>>>> RealMatrix Qcoeff = coeffQRd.getQ();
>>>>>> RealMatrix R = X.multiply(Qcoeff);
>>>>>>
>>>>>> final int nvars = X.getColumnDimension();
>>>>>> final int nobs = X.getRowDimension();
>>>>>> final int ncons = coeff.length;
>>>>>>
>>>>>> RealMatrix R2 = R.getSubMatrix(
>>>>>> 0, nobs - 1, ncons, nvars - 1);
>>>>>>
>>>>>> RealMatrix R1 = R.getSubMatrix(
>>>>>> 0, nobs - 1, 0, ncons - 1);
>>>>>>
>>>>>> RealVector gamma = rhsVec.copy();
>>>>>>
>>>>>> RealMatrix coeffR = coeffQRd.getR().getSubMatrix(
>>>>>> 0, ncons - 1, 0, ncons - 1);
>>>>>>
>>>>>> MatrixUtils.solveLowerTriangularSystem(coeffR.transpose(),
>>>> gamma);
>>>>>> RealVector gammPrime = Y.subtract(R1.operate(gamma));
>>>>>>
>>>>>> QRDecomposition qr2 = new QRDecompositionImpl(R2);
>>>>>>
>>>>>> RealVector constrainedSolution =
>>>>>> (qr2.getSolver().solve(gammPrime));
>>>>>>
>>>>>> RealVector stackedVector =
>>>>>> new ArrayRealVector(
>>>>>> gamma.toArray(),
>>>>>> constrainedSolution.toArray());
>>>>>>
>>>>>> stackedVector = Qcoeff.operate(stackedVector);
>>>>>>
>>>>>> return stackedVector.toArray();
>>>>>> }
>>>>>>
>>>>>> This approach is based on Dongarra et al:
>>>>>>
>>>>>> LAPACK Working Note
>>>>>> Generalized QR Factorization and its Applications
>>>>>> Work in Progress
>>>>>> E. Anderson, Z. Bai and J. Dongarra
>>>>>> December 9, 1991
>>>>>> August 9, 1994
>>>>>>
>>>>>> There is nothing terrible about this approach, the coding is not
>>>> finished
>>>>>> and tidy, but its a work in progress.
>>>>>>
>>>>>> I am also aware of second approach. I do not have a cite for it, I
>>> think
>>>> I
>>>>>> may have derived it myself, but it would not surprise me if it is in
>>>> some
>>>>>> textbook somewhere... That second approach takes the QR
>> decomposition
>>> of
>>>> the
>>>>>> coefficient matrix and calculates adjustment matrices for the design
>>>> matrix
>>>>>> and dependent vector. The problem is that I need to reorganize the
>>>> design
>>>>>> matrix by the pivots of the QR decomposition. Once I have the
>>> adjustment
>>>>>> matrices, everything should proceed as in the case of an
>> unconstrained
>>>>>> estimation. I like the idea that if we transform the data,
>> everything
>>>> works
>>>>>> the same way.
>>>>>>
>>>>>> Since then the ConstrainedOLSMultipleRegression class looks like:
>>>>>> public ConstrainedOLSMultipleRegression extends
>> OLSMultipleRegression{
>>>>>> }
>>>>>>
>>>>>>
>>>>>> As for the fact that the QRDecompositionImpl reflects its interface.
>>> We
>>>>>> should probably add the functions:
>>>>>> public int[] getPivots();
>>>>>> public boolean isPivotting();
>>>>>>
>>>>>> to the interface. As Christopher pointed out, if the current
>>>> decomposition
>>>>>> is non pivoting, its pivot record is the canonical one,
>>> {0,1,2,...,n-1}.
>>>>>> -Greg
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>
>>>> ---------------------------------------------------------------------
>>>> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
>>>> For additional commands, e-mail: dev-help@commons.apache.org
>>>>
>>>>
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] Pivoting in QR decomposition
Posted by Ted Dunning <te...@gmail.com>.
Also, if the QRDecomposition interface *is* extended with getP, it is the
work of a moment to change it to an abstract class instead of an interface
and provide a default method for getP that throws
UnsupportedOperationException. Since there are very unlikely to be any user
extensions of QRDecomposition in the wild, we might as well do both changes
at once. For that matter, getP can also easily return an identity
permutation by default.
Much like an immutable set or list implements the standard remove method,
but complains when it is called, it is quite reasonable for QRDecomposition
to not be quite the lowest common denominator, but something more like a
mid-line that contains the things that you might expect a QR decomposition
to do. Having a few (very few) methods of this sort is preferable to having
an elaborate class structure of alternative interfaces.
On Sat, Jul 23, 2011 at 12:45 PM, Chris Nix <ch...@gmail.com> wrote:
> We can do this with the current Householder reflection implementation,
> except instead of just obtaining reflections from columns in sequence
> across
> the input matrix, we select the column with the greatest L2-norm at each
> iteration. The resulting permutation matrix is thus built and can be
> returned later with a getP() method. A pleasing by-product is that the
> resulting R matrix is 'rank-revealing', allowing for a quicker getRank()
> than currently exists in the SingularValueDecompositionImpl class.
>
> Does it make sense to extend the current QRDecomposition interface to one
> for rank-revealing QR decompositions that have the existing methods, plus a
> getP() and getRank()? The implementing class could extend the current
> QRDecompositionImpl class to save reproducing code, at the cost of opening
> up some private variables and methods, and the solver.
>
> I'll open a JIRA issue.
>
> Chris.
>
> On 23 July 2011 18:44, Greg Sterijevski <gs...@gmail.com> wrote:
>
> > Chris, you had an algorithm in mind? -Greg
> >
> > On Sat, Jul 23, 2011 at 11:29 AM, Phil Steitz <ph...@gmail.com>
> > wrote:
> >
> > > On 7/22/11 11:40 AM, Greg Sterijevski wrote:
> > > > Sorry,
> > > >
> > > > public ConstrainedOLSMultipleRegression extends
> OLSMultipleRegression{}
> > > >
> > > > should read:
> > > >
> > > > public ConstrainedOLSMultipleRegression extends
> OLSMultipleRegression{
> > > > @Override
> > > > public void newSampleData(double[] data, double[][] coeff,
> double[]
> > > rhs,
> > > > int nob, int nvars) {
> > > > adjustData( data, coeff, rhs);
> > > > super.newSampleData(data, nobs, nvars);
> > > > qr = new QRDecompositionImpl(X);
> > > > }
> > > >
> > > >> }
> > > > The data would be transformed on the way in, and everything else
> would
> > > > remain the same...
> > >
> > > Got it. Sounds good. Patch away...
> > >
> > > Couple of things to keep in mind:
> > >
> > > 0) We may want to dispense with the QRDecomposition interface
> > > altogther. If we keep it, we should trim it down to what is common
> > > and meaningfully implemented in all impls. So both the Householder
> > > and permutation getters are dropped. If you need a pivoting impl,
> > > go a head and code it up and we can reassess the interface.
> > >
> > > 1) We should be aiming to standardize the regression API. Lets pick
> > > up the other thread on regression refactoring. Before hacking too
> > > much more on OLS, lets refine and retrofit the new regression API on
> > > this class.
> > >
> > > Phil
> > > >
> > > >
> > > > On Fri, Jul 22, 2011 at 1:37 PM, Greg Sterijevski <
> > > gsterijevski@gmail.com>wrote:
> > > >
> > > >> On the need for pivoting:
> > > >>
> > > >> Here is my first approach for changing OLSMultipleRegression to do
> > > >> constrained estimation:
> > > >>
> > > >> public double[] calculateBeta(double[][] coeff, double[] rhs) {
> > > >> if (rhs.length != coeff.length) {
> > > >> throw new IllegalArgumentException("");
> > > >> }
> > > >> for (double[] rest : coeff) {
> > > >> if (rest.length != this.X.getColumnDimension()) {
> > > >> throw new IllegalArgumentException("");
> > > >> }
> > > >> }
> > > >> RealMatrix Coeff = new Array2DRowRealMatrix(coeff, false);
> > > >> RealVector rhsVec = new ArrayRealVector(rhs);
> > > >> QRDecomposition coeffQRd = new
> > > >> QRDecompositionImpl(Coeff.transpose());
> > > >> RealMatrix Qcoeff = coeffQRd.getQ();
> > > >> RealMatrix R = X.multiply(Qcoeff);
> > > >>
> > > >> final int nvars = X.getColumnDimension();
> > > >> final int nobs = X.getRowDimension();
> > > >> final int ncons = coeff.length;
> > > >>
> > > >> RealMatrix R2 = R.getSubMatrix(
> > > >> 0, nobs - 1, ncons, nvars - 1);
> > > >>
> > > >> RealMatrix R1 = R.getSubMatrix(
> > > >> 0, nobs - 1, 0, ncons - 1);
> > > >>
> > > >> RealVector gamma = rhsVec.copy();
> > > >>
> > > >> RealMatrix coeffR = coeffQRd.getR().getSubMatrix(
> > > >> 0, ncons - 1, 0, ncons - 1);
> > > >>
> > > >> MatrixUtils.solveLowerTriangularSystem(coeffR.transpose(),
> > > gamma);
> > > >>
> > > >> RealVector gammPrime = Y.subtract(R1.operate(gamma));
> > > >>
> > > >> QRDecomposition qr2 = new QRDecompositionImpl(R2);
> > > >>
> > > >> RealVector constrainedSolution =
> > > >> (qr2.getSolver().solve(gammPrime));
> > > >>
> > > >> RealVector stackedVector =
> > > >> new ArrayRealVector(
> > > >> gamma.toArray(),
> > > >> constrainedSolution.toArray());
> > > >>
> > > >> stackedVector = Qcoeff.operate(stackedVector);
> > > >>
> > > >> return stackedVector.toArray();
> > > >> }
> > > >>
> > > >> This approach is based on Dongarra et al:
> > > >>
> > > >> LAPACK Working Note
> > > >> Generalized QR Factorization and its Applications
> > > >> Work in Progress
> > > >> E. Anderson, Z. Bai and J. Dongarra
> > > >> December 9, 1991
> > > >> August 9, 1994
> > > >>
> > > >> There is nothing terrible about this approach, the coding is not
> > > finished
> > > >> and tidy, but its a work in progress.
> > > >>
> > > >> I am also aware of second approach. I do not have a cite for it, I
> > think
> > > I
> > > >> may have derived it myself, but it would not surprise me if it is in
> > > some
> > > >> textbook somewhere... That second approach takes the QR
> decomposition
> > of
> > > the
> > > >> coefficient matrix and calculates adjustment matrices for the design
> > > matrix
> > > >> and dependent vector. The problem is that I need to reorganize the
> > > design
> > > >> matrix by the pivots of the QR decomposition. Once I have the
> > adjustment
> > > >> matrices, everything should proceed as in the case of an
> unconstrained
> > > >> estimation. I like the idea that if we transform the data,
> everything
> > > works
> > > >> the same way.
> > > >>
> > > >> Since then the ConstrainedOLSMultipleRegression class looks like:
> > > >> public ConstrainedOLSMultipleRegression extends
> OLSMultipleRegression{
> > > >>
> > > >> }
> > > >>
> > > >>
> > > >> As for the fact that the QRDecompositionImpl reflects its interface.
> > We
> > > >> should probably add the functions:
> > > >> public int[] getPivots();
> > > >> public boolean isPivotting();
> > > >>
> > > >> to the interface. As Christopher pointed out, if the current
> > > decomposition
> > > >> is non pivoting, its pivot record is the canonical one,
> > {0,1,2,...,n-1}.
> > > >>
> > > >> -Greg
> > > >>
> > > >>
> > > >>
> > > >>
> > > >>
> > >
> > >
> > > ---------------------------------------------------------------------
> > > To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> > > For additional commands, e-mail: dev-help@commons.apache.org
> > >
> > >
> >
>
Re: [math] Pivoting in QR decomposition
Posted by Chris Nix <ch...@gmail.com>.
We can do this with the current Householder reflection implementation,
except instead of just obtaining reflections from columns in sequence across
the input matrix, we select the column with the greatest L2-norm at each
iteration. The resulting permutation matrix is thus built and can be
returned later with a getP() method. A pleasing by-product is that the
resulting R matrix is 'rank-revealing', allowing for a quicker getRank()
than currently exists in the SingularValueDecompositionImpl class.
Does it make sense to extend the current QRDecomposition interface to one
for rank-revealing QR decompositions that have the existing methods, plus a
getP() and getRank()? The implementing class could extend the current
QRDecompositionImpl class to save reproducing code, at the cost of opening
up some private variables and methods, and the solver.
I'll open a JIRA issue.
Chris.
On 23 July 2011 18:44, Greg Sterijevski <gs...@gmail.com> wrote:
> Chris, you had an algorithm in mind? -Greg
>
> On Sat, Jul 23, 2011 at 11:29 AM, Phil Steitz <ph...@gmail.com>
> wrote:
>
> > On 7/22/11 11:40 AM, Greg Sterijevski wrote:
> > > Sorry,
> > >
> > > public ConstrainedOLSMultipleRegression extends OLSMultipleRegression{}
> > >
> > > should read:
> > >
> > > public ConstrainedOLSMultipleRegression extends OLSMultipleRegression{
> > > @Override
> > > public void newSampleData(double[] data, double[][] coeff, double[]
> > rhs,
> > > int nob, int nvars) {
> > > adjustData( data, coeff, rhs);
> > > super.newSampleData(data, nobs, nvars);
> > > qr = new QRDecompositionImpl(X);
> > > }
> > >
> > >> }
> > > The data would be transformed on the way in, and everything else would
> > > remain the same...
> >
> > Got it. Sounds good. Patch away...
> >
> > Couple of things to keep in mind:
> >
> > 0) We may want to dispense with the QRDecomposition interface
> > altogther. If we keep it, we should trim it down to what is common
> > and meaningfully implemented in all impls. So both the Householder
> > and permutation getters are dropped. If you need a pivoting impl,
> > go a head and code it up and we can reassess the interface.
> >
> > 1) We should be aiming to standardize the regression API. Lets pick
> > up the other thread on regression refactoring. Before hacking too
> > much more on OLS, lets refine and retrofit the new regression API on
> > this class.
> >
> > Phil
> > >
> > >
> > > On Fri, Jul 22, 2011 at 1:37 PM, Greg Sterijevski <
> > gsterijevski@gmail.com>wrote:
> > >
> > >> On the need for pivoting:
> > >>
> > >> Here is my first approach for changing OLSMultipleRegression to do
> > >> constrained estimation:
> > >>
> > >> public double[] calculateBeta(double[][] coeff, double[] rhs) {
> > >> if (rhs.length != coeff.length) {
> > >> throw new IllegalArgumentException("");
> > >> }
> > >> for (double[] rest : coeff) {
> > >> if (rest.length != this.X.getColumnDimension()) {
> > >> throw new IllegalArgumentException("");
> > >> }
> > >> }
> > >> RealMatrix Coeff = new Array2DRowRealMatrix(coeff, false);
> > >> RealVector rhsVec = new ArrayRealVector(rhs);
> > >> QRDecomposition coeffQRd = new
> > >> QRDecompositionImpl(Coeff.transpose());
> > >> RealMatrix Qcoeff = coeffQRd.getQ();
> > >> RealMatrix R = X.multiply(Qcoeff);
> > >>
> > >> final int nvars = X.getColumnDimension();
> > >> final int nobs = X.getRowDimension();
> > >> final int ncons = coeff.length;
> > >>
> > >> RealMatrix R2 = R.getSubMatrix(
> > >> 0, nobs - 1, ncons, nvars - 1);
> > >>
> > >> RealMatrix R1 = R.getSubMatrix(
> > >> 0, nobs - 1, 0, ncons - 1);
> > >>
> > >> RealVector gamma = rhsVec.copy();
> > >>
> > >> RealMatrix coeffR = coeffQRd.getR().getSubMatrix(
> > >> 0, ncons - 1, 0, ncons - 1);
> > >>
> > >> MatrixUtils.solveLowerTriangularSystem(coeffR.transpose(),
> > gamma);
> > >>
> > >> RealVector gammPrime = Y.subtract(R1.operate(gamma));
> > >>
> > >> QRDecomposition qr2 = new QRDecompositionImpl(R2);
> > >>
> > >> RealVector constrainedSolution =
> > >> (qr2.getSolver().solve(gammPrime));
> > >>
> > >> RealVector stackedVector =
> > >> new ArrayRealVector(
> > >> gamma.toArray(),
> > >> constrainedSolution.toArray());
> > >>
> > >> stackedVector = Qcoeff.operate(stackedVector);
> > >>
> > >> return stackedVector.toArray();
> > >> }
> > >>
> > >> This approach is based on Dongarra et al:
> > >>
> > >> LAPACK Working Note
> > >> Generalized QR Factorization and its Applications
> > >> Work in Progress
> > >> E. Anderson, Z. Bai and J. Dongarra
> > >> December 9, 1991
> > >> August 9, 1994
> > >>
> > >> There is nothing terrible about this approach, the coding is not
> > finished
> > >> and tidy, but its a work in progress.
> > >>
> > >> I am also aware of second approach. I do not have a cite for it, I
> think
> > I
> > >> may have derived it myself, but it would not surprise me if it is in
> > some
> > >> textbook somewhere... That second approach takes the QR decomposition
> of
> > the
> > >> coefficient matrix and calculates adjustment matrices for the design
> > matrix
> > >> and dependent vector. The problem is that I need to reorganize the
> > design
> > >> matrix by the pivots of the QR decomposition. Once I have the
> adjustment
> > >> matrices, everything should proceed as in the case of an unconstrained
> > >> estimation. I like the idea that if we transform the data, everything
> > works
> > >> the same way.
> > >>
> > >> Since then the ConstrainedOLSMultipleRegression class looks like:
> > >> public ConstrainedOLSMultipleRegression extends OLSMultipleRegression{
> > >>
> > >> }
> > >>
> > >>
> > >> As for the fact that the QRDecompositionImpl reflects its interface.
> We
> > >> should probably add the functions:
> > >> public int[] getPivots();
> > >> public boolean isPivotting();
> > >>
> > >> to the interface. As Christopher pointed out, if the current
> > decomposition
> > >> is non pivoting, its pivot record is the canonical one,
> {0,1,2,...,n-1}.
> > >>
> > >> -Greg
> > >>
> > >>
> > >>
> > >>
> > >>
> >
> >
> > ---------------------------------------------------------------------
> > To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> > For additional commands, e-mail: dev-help@commons.apache.org
> >
> >
>
Re: [math] Pivoting in QR decomposition
Posted by Greg Sterijevski <gs...@gmail.com>.
Chris, you had an algorithm in mind? -Greg
On Sat, Jul 23, 2011 at 11:29 AM, Phil Steitz <ph...@gmail.com> wrote:
> On 7/22/11 11:40 AM, Greg Sterijevski wrote:
> > Sorry,
> >
> > public ConstrainedOLSMultipleRegression extends OLSMultipleRegression{}
> >
> > should read:
> >
> > public ConstrainedOLSMultipleRegression extends OLSMultipleRegression{
> > @Override
> > public void newSampleData(double[] data, double[][] coeff, double[]
> rhs,
> > int nob, int nvars) {
> > adjustData( data, coeff, rhs);
> > super.newSampleData(data, nobs, nvars);
> > qr = new QRDecompositionImpl(X);
> > }
> >
> >> }
> > The data would be transformed on the way in, and everything else would
> > remain the same...
>
> Got it. Sounds good. Patch away...
>
> Couple of things to keep in mind:
>
> 0) We may want to dispense with the QRDecomposition interface
> altogther. If we keep it, we should trim it down to what is common
> and meaningfully implemented in all impls. So both the Householder
> and permutation getters are dropped. If you need a pivoting impl,
> go a head and code it up and we can reassess the interface.
>
> 1) We should be aiming to standardize the regression API. Lets pick
> up the other thread on regression refactoring. Before hacking too
> much more on OLS, lets refine and retrofit the new regression API on
> this class.
>
> Phil
> >
> >
> > On Fri, Jul 22, 2011 at 1:37 PM, Greg Sterijevski <
> gsterijevski@gmail.com>wrote:
> >
> >> On the need for pivoting:
> >>
> >> Here is my first approach for changing OLSMultipleRegression to do
> >> constrained estimation:
> >>
> >> public double[] calculateBeta(double[][] coeff, double[] rhs) {
> >> if (rhs.length != coeff.length) {
> >> throw new IllegalArgumentException("");
> >> }
> >> for (double[] rest : coeff) {
> >> if (rest.length != this.X.getColumnDimension()) {
> >> throw new IllegalArgumentException("");
> >> }
> >> }
> >> RealMatrix Coeff = new Array2DRowRealMatrix(coeff, false);
> >> RealVector rhsVec = new ArrayRealVector(rhs);
> >> QRDecomposition coeffQRd = new
> >> QRDecompositionImpl(Coeff.transpose());
> >> RealMatrix Qcoeff = coeffQRd.getQ();
> >> RealMatrix R = X.multiply(Qcoeff);
> >>
> >> final int nvars = X.getColumnDimension();
> >> final int nobs = X.getRowDimension();
> >> final int ncons = coeff.length;
> >>
> >> RealMatrix R2 = R.getSubMatrix(
> >> 0, nobs - 1, ncons, nvars - 1);
> >>
> >> RealMatrix R1 = R.getSubMatrix(
> >> 0, nobs - 1, 0, ncons - 1);
> >>
> >> RealVector gamma = rhsVec.copy();
> >>
> >> RealMatrix coeffR = coeffQRd.getR().getSubMatrix(
> >> 0, ncons - 1, 0, ncons - 1);
> >>
> >> MatrixUtils.solveLowerTriangularSystem(coeffR.transpose(),
> gamma);
> >>
> >> RealVector gammPrime = Y.subtract(R1.operate(gamma));
> >>
> >> QRDecomposition qr2 = new QRDecompositionImpl(R2);
> >>
> >> RealVector constrainedSolution =
> >> (qr2.getSolver().solve(gammPrime));
> >>
> >> RealVector stackedVector =
> >> new ArrayRealVector(
> >> gamma.toArray(),
> >> constrainedSolution.toArray());
> >>
> >> stackedVector = Qcoeff.operate(stackedVector);
> >>
> >> return stackedVector.toArray();
> >> }
> >>
> >> This approach is based on Dongarra et al:
> >>
> >> LAPACK Working Note
> >> Generalized QR Factorization and its Applications
> >> Work in Progress
> >> E. Anderson, Z. Bai and J. Dongarra
> >> December 9, 1991
> >> August 9, 1994
> >>
> >> There is nothing terrible about this approach, the coding is not
> finished
> >> and tidy, but its a work in progress.
> >>
> >> I am also aware of second approach. I do not have a cite for it, I think
> I
> >> may have derived it myself, but it would not surprise me if it is in
> some
> >> textbook somewhere... That second approach takes the QR decomposition of
> the
> >> coefficient matrix and calculates adjustment matrices for the design
> matrix
> >> and dependent vector. The problem is that I need to reorganize the
> design
> >> matrix by the pivots of the QR decomposition. Once I have the adjustment
> >> matrices, everything should proceed as in the case of an unconstrained
> >> estimation. I like the idea that if we transform the data, everything
> works
> >> the same way.
> >>
> >> Since then the ConstrainedOLSMultipleRegression class looks like:
> >> public ConstrainedOLSMultipleRegression extends OLSMultipleRegression{
> >>
> >> }
> >>
> >>
> >> As for the fact that the QRDecompositionImpl reflects its interface. We
> >> should probably add the functions:
> >> public int[] getPivots();
> >> public boolean isPivotting();
> >>
> >> to the interface. As Christopher pointed out, if the current
> decomposition
> >> is non pivoting, its pivot record is the canonical one, {0,1,2,...,n-1}.
> >>
> >> -Greg
> >>
> >>
> >>
> >>
> >>
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
>
Re: [math] Pivoting in QR decomposition
Posted by Greg Sterijevski <gs...@gmail.com>.
I agree about refactoring, was just experimenting to see what sorts of
changes would need to be made in order to support restrictions. In the next
week or so, I will introduce yet another regression algorithm, the so-called
"sweep" algorithm. Once we have that, we will cover both normal equation
based estimation techniques and build the inverse as you techniques. Then we
will have something concrete to finalize our refactoring discussions.
-Greg
On Sat, Jul 23, 2011 at 11:29 AM, Phil Steitz <ph...@gmail.com> wrote:
> On 7/22/11 11:40 AM, Greg Sterijevski wrote:
> > Sorry,
> >
> > public ConstrainedOLSMultipleRegression extends OLSMultipleRegression{}
> >
> > should read:
> >
> > public ConstrainedOLSMultipleRegression extends OLSMultipleRegression{
> > @Override
> > public void newSampleData(double[] data, double[][] coeff, double[]
> rhs,
> > int nob, int nvars) {
> > adjustData( data, coeff, rhs);
> > super.newSampleData(data, nobs, nvars);
> > qr = new QRDecompositionImpl(X);
> > }
> >
> >> }
> > The data would be transformed on the way in, and everything else would
> > remain the same...
>
> Got it. Sounds good. Patch away...
>
> Couple of things to keep in mind:
>
> 0) We may want to dispense with the QRDecomposition interface
> altogther. If we keep it, we should trim it down to what is common
> and meaningfully implemented in all impls. So both the Householder
> and permutation getters are dropped. If you need a pivoting impl,
> go a head and code it up and we can reassess the interface.
>
> 1) We should be aiming to standardize the regression API. Lets pick
> up the other thread on regression refactoring. Before hacking too
> much more on OLS, lets refine and retrofit the new regression API on
> this class.
>
> Phil
> >
> >
> > On Fri, Jul 22, 2011 at 1:37 PM, Greg Sterijevski <
> gsterijevski@gmail.com>wrote:
> >
> >> On the need for pivoting:
> >>
> >> Here is my first approach for changing OLSMultipleRegression to do
> >> constrained estimation:
> >>
> >> public double[] calculateBeta(double[][] coeff, double[] rhs) {
> >> if (rhs.length != coeff.length) {
> >> throw new IllegalArgumentException("");
> >> }
> >> for (double[] rest : coeff) {
> >> if (rest.length != this.X.getColumnDimension()) {
> >> throw new IllegalArgumentException("");
> >> }
> >> }
> >> RealMatrix Coeff = new Array2DRowRealMatrix(coeff, false);
> >> RealVector rhsVec = new ArrayRealVector(rhs);
> >> QRDecomposition coeffQRd = new
> >> QRDecompositionImpl(Coeff.transpose());
> >> RealMatrix Qcoeff = coeffQRd.getQ();
> >> RealMatrix R = X.multiply(Qcoeff);
> >>
> >> final int nvars = X.getColumnDimension();
> >> final int nobs = X.getRowDimension();
> >> final int ncons = coeff.length;
> >>
> >> RealMatrix R2 = R.getSubMatrix(
> >> 0, nobs - 1, ncons, nvars - 1);
> >>
> >> RealMatrix R1 = R.getSubMatrix(
> >> 0, nobs - 1, 0, ncons - 1);
> >>
> >> RealVector gamma = rhsVec.copy();
> >>
> >> RealMatrix coeffR = coeffQRd.getR().getSubMatrix(
> >> 0, ncons - 1, 0, ncons - 1);
> >>
> >> MatrixUtils.solveLowerTriangularSystem(coeffR.transpose(),
> gamma);
> >>
> >> RealVector gammPrime = Y.subtract(R1.operate(gamma));
> >>
> >> QRDecomposition qr2 = new QRDecompositionImpl(R2);
> >>
> >> RealVector constrainedSolution =
> >> (qr2.getSolver().solve(gammPrime));
> >>
> >> RealVector stackedVector =
> >> new ArrayRealVector(
> >> gamma.toArray(),
> >> constrainedSolution.toArray());
> >>
> >> stackedVector = Qcoeff.operate(stackedVector);
> >>
> >> return stackedVector.toArray();
> >> }
> >>
> >> This approach is based on Dongarra et al:
> >>
> >> LAPACK Working Note
> >> Generalized QR Factorization and its Applications
> >> Work in Progress
> >> E. Anderson, Z. Bai and J. Dongarra
> >> December 9, 1991
> >> August 9, 1994
> >>
> >> There is nothing terrible about this approach, the coding is not
> finished
> >> and tidy, but its a work in progress.
> >>
> >> I am also aware of second approach. I do not have a cite for it, I think
> I
> >> may have derived it myself, but it would not surprise me if it is in
> some
> >> textbook somewhere... That second approach takes the QR decomposition of
> the
> >> coefficient matrix and calculates adjustment matrices for the design
> matrix
> >> and dependent vector. The problem is that I need to reorganize the
> design
> >> matrix by the pivots of the QR decomposition. Once I have the adjustment
> >> matrices, everything should proceed as in the case of an unconstrained
> >> estimation. I like the idea that if we transform the data, everything
> works
> >> the same way.
> >>
> >> Since then the ConstrainedOLSMultipleRegression class looks like:
> >> public ConstrainedOLSMultipleRegression extends OLSMultipleRegression{
> >>
> >> }
> >>
> >>
> >> As for the fact that the QRDecompositionImpl reflects its interface. We
> >> should probably add the functions:
> >> public int[] getPivots();
> >> public boolean isPivotting();
> >>
> >> to the interface. As Christopher pointed out, if the current
> decomposition
> >> is non pivoting, its pivot record is the canonical one, {0,1,2,...,n-1}.
> >>
> >> -Greg
> >>
> >>
> >>
> >>
> >>
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
>
Re: [math] Pivoting in QR decomposition
Posted by Phil Steitz <ph...@gmail.com>.
On 7/22/11 11:40 AM, Greg Sterijevski wrote:
> Sorry,
>
> public ConstrainedOLSMultipleRegression extends OLSMultipleRegression{}
>
> should read:
>
> public ConstrainedOLSMultipleRegression extends OLSMultipleRegression{
> @Override
> public void newSampleData(double[] data, double[][] coeff, double[] rhs,
> int nob, int nvars) {
> adjustData( data, coeff, rhs);
> super.newSampleData(data, nobs, nvars);
> qr = new QRDecompositionImpl(X);
> }
>
>> }
> The data would be transformed on the way in, and everything else would
> remain the same...
Got it. Sounds good. Patch away...
Couple of things to keep in mind:
0) We may want to dispense with the QRDecomposition interface
altogther. If we keep it, we should trim it down to what is common
and meaningfully implemented in all impls. So both the Householder
and permutation getters are dropped. If you need a pivoting impl,
go a head and code it up and we can reassess the interface.
1) We should be aiming to standardize the regression API. Lets pick
up the other thread on regression refactoring. Before hacking too
much more on OLS, lets refine and retrofit the new regression API on
this class.
Phil
>
>
> On Fri, Jul 22, 2011 at 1:37 PM, Greg Sterijevski <gs...@gmail.com>wrote:
>
>> On the need for pivoting:
>>
>> Here is my first approach for changing OLSMultipleRegression to do
>> constrained estimation:
>>
>> public double[] calculateBeta(double[][] coeff, double[] rhs) {
>> if (rhs.length != coeff.length) {
>> throw new IllegalArgumentException("");
>> }
>> for (double[] rest : coeff) {
>> if (rest.length != this.X.getColumnDimension()) {
>> throw new IllegalArgumentException("");
>> }
>> }
>> RealMatrix Coeff = new Array2DRowRealMatrix(coeff, false);
>> RealVector rhsVec = new ArrayRealVector(rhs);
>> QRDecomposition coeffQRd = new
>> QRDecompositionImpl(Coeff.transpose());
>> RealMatrix Qcoeff = coeffQRd.getQ();
>> RealMatrix R = X.multiply(Qcoeff);
>>
>> final int nvars = X.getColumnDimension();
>> final int nobs = X.getRowDimension();
>> final int ncons = coeff.length;
>>
>> RealMatrix R2 = R.getSubMatrix(
>> 0, nobs - 1, ncons, nvars - 1);
>>
>> RealMatrix R1 = R.getSubMatrix(
>> 0, nobs - 1, 0, ncons - 1);
>>
>> RealVector gamma = rhsVec.copy();
>>
>> RealMatrix coeffR = coeffQRd.getR().getSubMatrix(
>> 0, ncons - 1, 0, ncons - 1);
>>
>> MatrixUtils.solveLowerTriangularSystem(coeffR.transpose(), gamma);
>>
>> RealVector gammPrime = Y.subtract(R1.operate(gamma));
>>
>> QRDecomposition qr2 = new QRDecompositionImpl(R2);
>>
>> RealVector constrainedSolution =
>> (qr2.getSolver().solve(gammPrime));
>>
>> RealVector stackedVector =
>> new ArrayRealVector(
>> gamma.toArray(),
>> constrainedSolution.toArray());
>>
>> stackedVector = Qcoeff.operate(stackedVector);
>>
>> return stackedVector.toArray();
>> }
>>
>> This approach is based on Dongarra et al:
>>
>> LAPACK Working Note
>> Generalized QR Factorization and its Applications
>> Work in Progress
>> E. Anderson, Z. Bai and J. Dongarra
>> December 9, 1991
>> August 9, 1994
>>
>> There is nothing terrible about this approach, the coding is not finished
>> and tidy, but its a work in progress.
>>
>> I am also aware of second approach. I do not have a cite for it, I think I
>> may have derived it myself, but it would not surprise me if it is in some
>> textbook somewhere... That second approach takes the QR decomposition of the
>> coefficient matrix and calculates adjustment matrices for the design matrix
>> and dependent vector. The problem is that I need to reorganize the design
>> matrix by the pivots of the QR decomposition. Once I have the adjustment
>> matrices, everything should proceed as in the case of an unconstrained
>> estimation. I like the idea that if we transform the data, everything works
>> the same way.
>>
>> Since then the ConstrainedOLSMultipleRegression class looks like:
>> public ConstrainedOLSMultipleRegression extends OLSMultipleRegression{
>>
>> }
>>
>>
>> As for the fact that the QRDecompositionImpl reflects its interface. We
>> should probably add the functions:
>> public int[] getPivots();
>> public boolean isPivotting();
>>
>> to the interface. As Christopher pointed out, if the current decomposition
>> is non pivoting, its pivot record is the canonical one, {0,1,2,...,n-1}.
>>
>> -Greg
>>
>>
>>
>>
>>
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] Pivoting in QR decomposition
Posted by Greg Sterijevski <gs...@gmail.com>.
Sorry,
public ConstrainedOLSMultipleRegression extends OLSMultipleRegression{}
should read:
public ConstrainedOLSMultipleRegression extends OLSMultipleRegression{
@Override
public void newSampleData(double[] data, double[][] coeff, double[] rhs,
int nob, int nvars) {
adjustData( data, coeff, rhs);
super.newSampleData(data, nobs, nvars);
qr = new QRDecompositionImpl(X);
}
>
> }
The data would be transformed on the way in, and everything else would
remain the same...
On Fri, Jul 22, 2011 at 1:37 PM, Greg Sterijevski <gs...@gmail.com>wrote:
> On the need for pivoting:
>
> Here is my first approach for changing OLSMultipleRegression to do
> constrained estimation:
>
> public double[] calculateBeta(double[][] coeff, double[] rhs) {
> if (rhs.length != coeff.length) {
> throw new IllegalArgumentException("");
> }
> for (double[] rest : coeff) {
> if (rest.length != this.X.getColumnDimension()) {
> throw new IllegalArgumentException("");
> }
> }
> RealMatrix Coeff = new Array2DRowRealMatrix(coeff, false);
> RealVector rhsVec = new ArrayRealVector(rhs);
> QRDecomposition coeffQRd = new
> QRDecompositionImpl(Coeff.transpose());
> RealMatrix Qcoeff = coeffQRd.getQ();
> RealMatrix R = X.multiply(Qcoeff);
>
> final int nvars = X.getColumnDimension();
> final int nobs = X.getRowDimension();
> final int ncons = coeff.length;
>
> RealMatrix R2 = R.getSubMatrix(
> 0, nobs - 1, ncons, nvars - 1);
>
> RealMatrix R1 = R.getSubMatrix(
> 0, nobs - 1, 0, ncons - 1);
>
> RealVector gamma = rhsVec.copy();
>
> RealMatrix coeffR = coeffQRd.getR().getSubMatrix(
> 0, ncons - 1, 0, ncons - 1);
>
> MatrixUtils.solveLowerTriangularSystem(coeffR.transpose(), gamma);
>
> RealVector gammPrime = Y.subtract(R1.operate(gamma));
>
> QRDecomposition qr2 = new QRDecompositionImpl(R2);
>
> RealVector constrainedSolution =
> (qr2.getSolver().solve(gammPrime));
>
> RealVector stackedVector =
> new ArrayRealVector(
> gamma.toArray(),
> constrainedSolution.toArray());
>
> stackedVector = Qcoeff.operate(stackedVector);
>
> return stackedVector.toArray();
> }
>
> This approach is based on Dongarra et al:
>
> LAPACK Working Note
> Generalized QR Factorization and its Applications
> Work in Progress
> E. Anderson, Z. Bai and J. Dongarra
> December 9, 1991
> August 9, 1994
>
> There is nothing terrible about this approach, the coding is not finished
> and tidy, but its a work in progress.
>
> I am also aware of second approach. I do not have a cite for it, I think I
> may have derived it myself, but it would not surprise me if it is in some
> textbook somewhere... That second approach takes the QR decomposition of the
> coefficient matrix and calculates adjustment matrices for the design matrix
> and dependent vector. The problem is that I need to reorganize the design
> matrix by the pivots of the QR decomposition. Once I have the adjustment
> matrices, everything should proceed as in the case of an unconstrained
> estimation. I like the idea that if we transform the data, everything works
> the same way.
>
> Since then the ConstrainedOLSMultipleRegression class looks like:
> public ConstrainedOLSMultipleRegression extends OLSMultipleRegression{
>
> }
>
>
> As for the fact that the QRDecompositionImpl reflects its interface. We
> should probably add the functions:
> public int[] getPivots();
> public boolean isPivotting();
>
> to the interface. As Christopher pointed out, if the current decomposition
> is non pivoting, its pivot record is the canonical one, {0,1,2,...,n-1}.
>
> -Greg
>
>
>
>
>
Re: [math] Pivoting in QR decomposition
Posted by Greg Sterijevski <gs...@gmail.com>.
On the need for pivoting:
Here is my first approach for changing OLSMultipleRegression to do
constrained estimation:
public double[] calculateBeta(double[][] coeff, double[] rhs) {
if (rhs.length != coeff.length) {
throw new IllegalArgumentException("");
}
for (double[] rest : coeff) {
if (rest.length != this.X.getColumnDimension()) {
throw new IllegalArgumentException("");
}
}
RealMatrix Coeff = new Array2DRowRealMatrix(coeff, false);
RealVector rhsVec = new ArrayRealVector(rhs);
QRDecomposition coeffQRd = new
QRDecompositionImpl(Coeff.transpose());
RealMatrix Qcoeff = coeffQRd.getQ();
RealMatrix R = X.multiply(Qcoeff);
final int nvars = X.getColumnDimension();
final int nobs = X.getRowDimension();
final int ncons = coeff.length;
RealMatrix R2 = R.getSubMatrix(
0, nobs - 1, ncons, nvars - 1);
RealMatrix R1 = R.getSubMatrix(
0, nobs - 1, 0, ncons - 1);
RealVector gamma = rhsVec.copy();
RealMatrix coeffR = coeffQRd.getR().getSubMatrix(
0, ncons - 1, 0, ncons - 1);
MatrixUtils.solveLowerTriangularSystem(coeffR.transpose(), gamma);
RealVector gammPrime = Y.subtract(R1.operate(gamma));
QRDecomposition qr2 = new QRDecompositionImpl(R2);
RealVector constrainedSolution = (qr2.getSolver().solve(gammPrime));
RealVector stackedVector =
new ArrayRealVector(
gamma.toArray(),
constrainedSolution.toArray());
stackedVector = Qcoeff.operate(stackedVector);
return stackedVector.toArray();
}
This approach is based on Dongarra et al:
LAPACK Working Note
Generalized QR Factorization and its Applications
Work in Progress
E. Anderson, Z. Bai and J. Dongarra
December 9, 1991
August 9, 1994
There is nothing terrible about this approach, the coding is not finished
and tidy, but its a work in progress.
I am also aware of second approach. I do not have a cite for it, I think I
may have derived it myself, but it would not surprise me if it is in some
textbook somewhere... That second approach takes the QR decomposition of the
coefficient matrix and calculates adjustment matrices for the design matrix
and dependent vector. The problem is that I need to reorganize the design
matrix by the pivots of the QR decomposition. Once I have the adjustment
matrices, everything should proceed as in the case of an unconstrained
estimation. I like the idea that if we transform the data, everything works
the same way.
Since then the ConstrainedOLSMultipleRegression class looks like:
public ConstrainedOLSMultipleRegression extends OLSMultipleRegression{
}
As for the fact that the QRDecompositionImpl reflects its interface. We
should probably add the functions:
public int[] getPivots();
public boolean isPivotting();
to the interface. As Christopher pointed out, if the current decomposition
is non pivoting, its pivot record is the canonical one, {0,1,2,...,n-1}.
-Greg
Re: [math] Pivoting in QR decomposition
Posted by Phil Steitz <ph...@gmail.com>.
On 7/21/11 6:56 PM, Greg Sterijevski wrote:
> If a pivoting QR decomp was to be included, I see it as existing in addition
> to the current nonpivoting version. I don't see this as an addition to the
> current QRDecompositionImpl. We would not be exposing more of the inner
> workings of the QRDecomposition, just a different view?
OK, if needed we can certainly add it. If and when we do that we
will have to deal with the fact that the current QRDecomposition
interface reflects QRDecompositionImpl.
Phil
>
> -Greg
>
> On Thu, Jul 21, 2011 at 6:39 PM, Ted Dunning <te...@gmail.com> wrote:
>
>> Also, pivoting is a good thing in a general QR decomposition since it
>> extends the applicability of the code to nearly rank deficient inputs.
>>
>> On Thu, Jul 21, 2011 at 3:17 PM, Greg Sterijevski <gsterijevski@gmail.com
>>> wrote:
>>> Sorry, did not mean to touch off a debate. There are two ways (that I
>> know
>>> of) to project linear restrictions through the data when running a
>>> regression using the QR decomposition. I need to review my notes, but I
>>> distinctly remember needing to get the pivots so that I know which
>> columns
>>> are 'contaminated' by the restriction... Also, in dealing with
>>> restrictions,
>>> I will need the pivots to clean up redundant restrictions. Doing a QR
>>> reduction on an arbitrary coefficient matrix will give me a reduced set
>> of
>>> restrictions and the restrictions which are mapped to the nullity (the
>>> redundancies).
>>>
>>> In any event, let me finish my experimentation and then I can more
>> cleanly
>>> enunciate my needs.
>>>
>>> On Thu, Jul 21, 2011 at 10:33 AM, Phil Steitz <ph...@gmail.com>
>>> wrote:
>>>
>>>> On 7/21/11 8:07 AM, Ted Dunning wrote:
>>>>> This is upside down.
>>>>>
>>>>> Opening the JIRA and putting the patch up is the best way to
>> determine
>>> if
>>>>> this is useful.
>>>>>
>>>>> On Thu, Jul 21, 2011 at 3:17 AM, Chris Nix <ch...@gmail.com>
>>> wrote:
>>>>>> I have just written a small patch that does QR Decomposition with
>>> column
>>>>>> pivoting and has a getP() method to get the resulting permutation
>>>> matrix.
>>>>>> If this is useful, I can open a JIRA issue for discussion.
>>>>>>
>>>> I think its best to understand more clearly what Greg's problem is
>>>> before taking time to prepare a patch. We have a working QR decomp
>>>> using Householder reflection. If it is useful to add another one,
>>>> we can talk about that; but I would like to understand more clearly
>>>> what, if any, problems the current implementation (and interface) has.
>>>>
>>>> Phil
>>>>
>>>> ---------------------------------------------------------------------
>>>> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
>>>> For additional commands, e-mail: dev-help@commons.apache.org
>>>>
>>>>
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] Pivoting in QR decomposition
Posted by Greg Sterijevski <gs...@gmail.com>.
If a pivoting QR decomp was to be included, I see it as existing in addition
to the current nonpivoting version. I don't see this as an addition to the
current QRDecompositionImpl. We would not be exposing more of the inner
workings of the QRDecomposition, just a different view?
-Greg
On Thu, Jul 21, 2011 at 6:39 PM, Ted Dunning <te...@gmail.com> wrote:
> Also, pivoting is a good thing in a general QR decomposition since it
> extends the applicability of the code to nearly rank deficient inputs.
>
> On Thu, Jul 21, 2011 at 3:17 PM, Greg Sterijevski <gsterijevski@gmail.com
> >wrote:
>
> > Sorry, did not mean to touch off a debate. There are two ways (that I
> know
> > of) to project linear restrictions through the data when running a
> > regression using the QR decomposition. I need to review my notes, but I
> > distinctly remember needing to get the pivots so that I know which
> columns
> > are 'contaminated' by the restriction... Also, in dealing with
> > restrictions,
> > I will need the pivots to clean up redundant restrictions. Doing a QR
> > reduction on an arbitrary coefficient matrix will give me a reduced set
> of
> > restrictions and the restrictions which are mapped to the nullity (the
> > redundancies).
> >
> > In any event, let me finish my experimentation and then I can more
> cleanly
> > enunciate my needs.
> >
> > On Thu, Jul 21, 2011 at 10:33 AM, Phil Steitz <ph...@gmail.com>
> > wrote:
> >
> > > On 7/21/11 8:07 AM, Ted Dunning wrote:
> > > > This is upside down.
> > > >
> > > > Opening the JIRA and putting the patch up is the best way to
> determine
> > if
> > > > this is useful.
> > > >
> > > > On Thu, Jul 21, 2011 at 3:17 AM, Chris Nix <ch...@gmail.com>
> > wrote:
> > > >
> > > >> I have just written a small patch that does QR Decomposition with
> > column
> > > >> pivoting and has a getP() method to get the resulting permutation
> > > matrix.
> > > >> If this is useful, I can open a JIRA issue for discussion.
> > > >>
> > > I think its best to understand more clearly what Greg's problem is
> > > before taking time to prepare a patch. We have a working QR decomp
> > > using Householder reflection. If it is useful to add another one,
> > > we can talk about that; but I would like to understand more clearly
> > > what, if any, problems the current implementation (and interface) has.
> > >
> > > Phil
> > >
> > > ---------------------------------------------------------------------
> > > To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> > > For additional commands, e-mail: dev-help@commons.apache.org
> > >
> > >
> >
>
Re: [math] Pivoting in QR decomposition
Posted by Ted Dunning <te...@gmail.com>.
Also, pivoting is a good thing in a general QR decomposition since it
extends the applicability of the code to nearly rank deficient inputs.
On Thu, Jul 21, 2011 at 3:17 PM, Greg Sterijevski <gs...@gmail.com>wrote:
> Sorry, did not mean to touch off a debate. There are two ways (that I know
> of) to project linear restrictions through the data when running a
> regression using the QR decomposition. I need to review my notes, but I
> distinctly remember needing to get the pivots so that I know which columns
> are 'contaminated' by the restriction... Also, in dealing with
> restrictions,
> I will need the pivots to clean up redundant restrictions. Doing a QR
> reduction on an arbitrary coefficient matrix will give me a reduced set of
> restrictions and the restrictions which are mapped to the nullity (the
> redundancies).
>
> In any event, let me finish my experimentation and then I can more cleanly
> enunciate my needs.
>
> On Thu, Jul 21, 2011 at 10:33 AM, Phil Steitz <ph...@gmail.com>
> wrote:
>
> > On 7/21/11 8:07 AM, Ted Dunning wrote:
> > > This is upside down.
> > >
> > > Opening the JIRA and putting the patch up is the best way to determine
> if
> > > this is useful.
> > >
> > > On Thu, Jul 21, 2011 at 3:17 AM, Chris Nix <ch...@gmail.com>
> wrote:
> > >
> > >> I have just written a small patch that does QR Decomposition with
> column
> > >> pivoting and has a getP() method to get the resulting permutation
> > matrix.
> > >> If this is useful, I can open a JIRA issue for discussion.
> > >>
> > I think its best to understand more clearly what Greg's problem is
> > before taking time to prepare a patch. We have a working QR decomp
> > using Householder reflection. If it is useful to add another one,
> > we can talk about that; but I would like to understand more clearly
> > what, if any, problems the current implementation (and interface) has.
> >
> > Phil
> >
> > ---------------------------------------------------------------------
> > To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> > For additional commands, e-mail: dev-help@commons.apache.org
> >
> >
>
Re: [math] Pivoting in QR decomposition
Posted by Greg Sterijevski <gs...@gmail.com>.
Sorry, did not mean to touch off a debate. There are two ways (that I know
of) to project linear restrictions through the data when running a
regression using the QR decomposition. I need to review my notes, but I
distinctly remember needing to get the pivots so that I know which columns
are 'contaminated' by the restriction... Also, in dealing with restrictions,
I will need the pivots to clean up redundant restrictions. Doing a QR
reduction on an arbitrary coefficient matrix will give me a reduced set of
restrictions and the restrictions which are mapped to the nullity (the
redundancies).
In any event, let me finish my experimentation and then I can more cleanly
enunciate my needs.
On Thu, Jul 21, 2011 at 10:33 AM, Phil Steitz <ph...@gmail.com> wrote:
> On 7/21/11 8:07 AM, Ted Dunning wrote:
> > This is upside down.
> >
> > Opening the JIRA and putting the patch up is the best way to determine if
> > this is useful.
> >
> > On Thu, Jul 21, 2011 at 3:17 AM, Chris Nix <ch...@gmail.com> wrote:
> >
> >> I have just written a small patch that does QR Decomposition with column
> >> pivoting and has a getP() method to get the resulting permutation
> matrix.
> >> If this is useful, I can open a JIRA issue for discussion.
> >>
> I think its best to understand more clearly what Greg's problem is
> before taking time to prepare a patch. We have a working QR decomp
> using Householder reflection. If it is useful to add another one,
> we can talk about that; but I would like to understand more clearly
> what, if any, problems the current implementation (and interface) has.
>
> Phil
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
>
Re: [math] Pivoting in QR decomposition
Posted by Phil Steitz <ph...@gmail.com>.
On 7/21/11 8:07 AM, Ted Dunning wrote:
> This is upside down.
>
> Opening the JIRA and putting the patch up is the best way to determine if
> this is useful.
>
> On Thu, Jul 21, 2011 at 3:17 AM, Chris Nix <ch...@gmail.com> wrote:
>
>> I have just written a small patch that does QR Decomposition with column
>> pivoting and has a getP() method to get the resulting permutation matrix.
>> If this is useful, I can open a JIRA issue for discussion.
>>
I think its best to understand more clearly what Greg's problem is
before taking time to prepare a patch. We have a working QR decomp
using Householder reflection. If it is useful to add another one,
we can talk about that; but I would like to understand more clearly
what, if any, problems the current implementation (and interface) has.
Phil
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org
Re: [math] Pivoting in QR decomposition
Posted by Ted Dunning <te...@gmail.com>.
This is upside down.
Opening the JIRA and putting the patch up is the best way to determine if
this is useful.
On Thu, Jul 21, 2011 at 3:17 AM, Chris Nix <ch...@gmail.com> wrote:
> I have just written a small patch that does QR Decomposition with column
> pivoting and has a getP() method to get the resulting permutation matrix.
> If this is useful, I can open a JIRA issue for discussion.
>
Re: [math] Pivoting in QR decomposition
Posted by Chris Nix <ch...@gmail.com>.
Looking at the current implementation, there is no column pivoting and so
the pivot matrix is the identity.
I have just written a small patch that does QR Decomposition with column
pivoting and has a getP() method to get the resulting permutation matrix.
If this is useful, I can open a JIRA issue for discussion.
Chris
On 21 July 2011 06:40, Phil Steitz <ph...@gmail.com> wrote:
> On 7/20/11 9:27 PM, Greg Sterijevski wrote:
> > Hello,
> >
> > In attempting to build constrained ols into the
> OLSMultipleLinearRegression
> > class, I rediscovered the fact that to keep the regression coefficients
> in
> > their canonical order, I need to know how the QR decomposition pivoted.
> In
> > other words, while I get the correct constrained parameters, their order
> is
> > changed. Is there any way to get this info out (pivot order)?
>
> Can you explain a little more exactly what you are trying to do?
> How exactly is the parameter vector getting permuted? The QRDecomp
> solver certainly should not do that. See the QRDecomp sources for
> what is exposed in the public API. We can talk about adding more,
> but we should be careful about exposing more implementation details.
> You can now get the Householder reflection vectors and obviously the
> decomp itself. As you can see from the sources, the current OLS
> impl just uses the solver.
>
> Phil
> >
> > Thanks,
> >
> > -Greg
> >
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
>
Re: [math] Pivoting in QR decomposition
Posted by Phil Steitz <ph...@gmail.com>.
On 7/20/11 9:27 PM, Greg Sterijevski wrote:
> Hello,
>
> In attempting to build constrained ols into the OLSMultipleLinearRegression
> class, I rediscovered the fact that to keep the regression coefficients in
> their canonical order, I need to know how the QR decomposition pivoted. In
> other words, while I get the correct constrained parameters, their order is
> changed. Is there any way to get this info out (pivot order)?
Can you explain a little more exactly what you are trying to do?
How exactly is the parameter vector getting permuted? The QRDecomp
solver certainly should not do that. See the QRDecomp sources for
what is exposed in the public API. We can talk about adding more,
but we should be careful about exposing more implementation details.
You can now get the Householder reflection vectors and obviously the
decomp itself. As you can see from the sources, the current OLS
impl just uses the solver.
Phil
>
> Thanks,
>
> -Greg
>
---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org