commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Ted Dunning <ted.dunn...@gmail.com>
Subject Re: [math] Pivoting in QR decomposition
Date Sat, 23 Jul 2011 20:37:35 GMT
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 <chris.nix@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 <gsterijevski@gmail.com> wrote:
>
> > Chris, you had an algorithm in mind? -Greg
> >
> > On Sat, Jul 23, 2011 at 11:29 AM, Phil Steitz <phil.steitz@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
> > >
> > >
> >
>

Mime
  • Unnamed multipart/alternative (inline, None, 0 bytes)
View raw message