I second +1. On Sat, Jul 23, 2011 at 4:06 PM, Phil Steitz <phil.ste...@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 <chris....@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 <gsterijev...@gmail.com> wrote: > >> > >>> Chris, you had an algorithm in mind? -Greg > >>> > >>> On Sat, Jul 23, 2011 at 11:29 AM, Phil Steitz <phil.ste...@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 < > >>>> gsterijev...@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-unsubscr...@commons.apache.org > >>>> For additional commands, e-mail: dev-h...@commons.apache.org > >>>> > >>>> > > > --------------------------------------------------------------------- > To unsubscribe, e-mail: dev-unsubscr...@commons.apache.org > For additional commands, e-mail: dev-h...@commons.apache.org > >