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
