@@ -92,6 +92,8 @@ public class QRDecomposition {
9292 */
9393 private DoubleMatrix2D Qcached = null ;
9494
95+ private DoubleMatrix2D effects ;
96+
9597 /**
9698 * @param A the matrix to decompose, pivoting will be used.
9799 */
@@ -227,39 +229,41 @@ public DoubleMatrix2D chol2inv() {
227229 }
228230
229231 /**
230- * Compute effects matrix Q'y (as in Rb = Q'y). Returns only
231- * the aspect associated with the parameters - the first <tt>rank</tt> elements.
232+ * Compute effects matrix Q'y (as in Rb = Q'y).
233+ *
234+ * <p>
235+ * "Tthe effects are the uncorrelated single-degree-of-freedom values obtained by projecting the data onto the
236+ * successive orthogonal subspaces generated by the QR decomposition during the fitting process. The first r (the
237+ * rank of the model) are associated with coefficients and the remainder span the space of residuals (but are not
238+ * associated with particular residuals)."
232239 *
233240 * @param y vector Missing values are ignored, otherwise assumed to be of the right size
234241 * @return vector of effects - these are the projections of y into Q column space
235242 */
236243 public DoubleMatrix1D effects ( DoubleMatrix1D y ) {
237244
238- Algebra solver = new Algebra ();
239- // Inefficient due to getting Q explicitly
240- return solver .mult ( solver .transpose ( this .getQ () ), MatrixUtil .removeMissing ( y ) ).viewPart ( 0 , this .rank );
241-
242- // Could do it this way, but QR.toArray() makes a copy...
243- // double[] qty = new double[y.size()];
244- // double[] junk = new double[y.size()];
245- // ubic.basecode.math.linalg.Dqrsl.dqrsl_j( QR.toArray(), QR.rows(), QR.columns(), qraux.toArray(), y.toArray(),
246- // junk, qty,
247- // junk, junk, junk, 1000 );
248- // return new DenseDoubleMatrix1D( qty );
245+ double [] qty = new double [y .size ()];
246+ double [] junk = new double [y .size ()];
247+ ubic .basecode .math .linalg .Dqrsl .dqrsl_j ( QR .toArray (), QR .rows (), QR .columns (), qraux .toArray (), MatrixUtil .removeMissing ( y ).toArray (),
248+ junk , qty ,
249+ junk , junk , junk , 1000 );
250+ return new DenseDoubleMatrix1D ( qty );
249251 }
250252
251253 /**
252- * Compute effects matrix Q'y (as in Rb = Q'y) Returns only the aspect associated with the parameters - the first
253- * <tt>rank</tt> elements.
254+ * Compute effects matrix Q'y (as in Rb = Q'y)
254255 *
255256 * @param y matrix of data, assumed to be of right size, missing values are not supported
256257 * @return matrix of effects - these are the projections of y's columns into Q column subspace associated with the
257258 * parameters,
258259 * so values are "effects" each basis vector on the data
259260 */
260261 public DoubleMatrix2D effects ( DoubleMatrix2D y ) {
261- Algebra solver = new Algebra ();
262- return solver .mult ( solver .transpose ( this .getQ () ), y );
262+ double [][] efa = new double [y .columns ()][y .rows ()];
263+ for ( int i = 0 ; i < y .columns (); i ++ ) {
264+ efa [i ] = effects ( y .viewColumn ( i ) ).toArray ();
265+ }
266+ return new DenseDoubleMatrix2D ( efa ).viewDice ();
263267 }
264268
265269 /**
@@ -387,11 +391,7 @@ public DoubleMatrix2D solve( DoubleMatrix2D y ) {
387391 throw new IllegalArgumentException ( "Matrix is rank deficient; try using pivoting" );
388392 }
389393
390- /*
391- * Direct computation of Q'y can be avoided by manipulation of the compact QR instead (see dqrsl).
392- *
393- */
394- DoubleMatrix2D qTy = effects ( y ); // will be modified by this call.
394+ DoubleMatrix2D qTy = effects ( y ); // FIXME we use this again later, but we recompute it. Try to cache it.
395395
396396 // Solve R*X = Y => X = RinvY; backsubstitution
397397 for ( int k1 = rank - 1 ; k1 >= 0 ; k1 -- ) {
@@ -421,7 +421,7 @@ public DoubleMatrix2D solve( DoubleMatrix2D y ) {
421421 */
422422 DoubleMatrix2D coeff = qTy .like ( p , y .columns () );
423423 coeff .assign ( Double .NaN );
424- for ( int i = 0 ; i < qTy . rows () ; i ++ ) {
424+ for ( int i = 0 ; i < this . rank ; i ++ ) {
425425 int piv = jpvt [i ]; // where the value should go.
426426 for ( int j = 0 ; j < qTy .columns (); j ++ ) {
427427 coeff .setQuick ( piv , j , qTy .getQuick ( i , j ) );
0 commit comments