commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Phil Steitz <phil.ste...@gmail.com>
Subject Re: [math] Pivoting in QR decomposition
Date Sat, 23 Jul 2011 21:06:51 GMT
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.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
>>>>
>>>>


---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org


Mime
View raw message