A bit of background to my request. I am writing interfaces to Eigen from R (http://www.R-project.org), a system for statistical computing, using Rcpp (CRAN.R-project.org/web/packages/Rcpp/), a C++ template library for the R internal storage types.

As a first example for such interfaces we usually create a method for least squares fits, including some of the statistical summaries. My current version of the code is

Code: Select all

extern "C" SEXP fastLm(SEXP ys, SEXP Xs) {
    using namespace Eigen;
    using namespace Rcpp;
    try {
   NumericMatrix X(Xs);
   NumericVector y(ys);
   size_t n = X.nrow(), p = X.ncol();
   if ((size_t)y.size() != n)
       throw std::invalid_argument("size mismatch");
   
   MatrixXd A = Map<MatrixXd>(X.begin(), n, p); // shares storage
   VectorXd b = Map<VectorXd>(y.begin(), n);
   
   ColPivHouseholderQR<MatrixXd> mqr(A);      // decompose the model matrix
   size_t               r = mqr.rank();
   ptrdiff_t           df = n - r;
   VectorXd          coef = mqr.solve(b);
   VectorXd        fitted = A * coef;
   VectorXd           res = b - fitted;
   double               s = std::sqrt(res.squaredNorm()/df);

   std::cout << "n = " << n << ", p = " << p << ", r = " << r << ", df = " << df << std::endl;
   MatrixXd  R = mqr.matrixQR().topLeftCorner(r, r); // need to create this explicitly
   std::cout << R << std::endl;
   MatrixXd  Rinv(TriangularView<MatrixXd, Upper>(R).solve(MatrixXd::Identity(r, r)));
   std::cout << Rinv << std::endl;

   VectorXd se = mqr.colsPermutation() * Rinv.rowwise().norm() * s; // standard errors
            // These parts can be replaced by wrap, once that is written
   NumericVector ccoef(coef.size()), rres(res.size()), sse(se.size()), ffitted(fitted.size());
   std::copy(coef.data(), coef.data() + coef.size(), ccoef.begin());
   std::copy(res.data(), res.data() + res.size(), rres.begin());
   std::copy(se.data(), se.data() + se.size(), sse.begin());
   std::copy(fitted.data(), fitted.data() + fitted.size(), ffitted.begin());
   IntegerVector perm(p);
   const int *pind = mqr.colsPermutation().indices().data();
   std::copy(pind, pind + p, perm.begin());

   return List::create(_["coefficients"]  = ccoef,
             _["rank"]          = r,
             _["df"]            = df,
             _["perm"]          = perm,
             _["stderr"]        = sse,
             _["s"]             = s,
             _["residuals"]     = rres,
             _["fitted.values"] = ffitted,
             _["Rinv"]          = NumericMatrix(r, r, Rinv.data())
       );
    } catch( std::exception &ex ) {
   forward_exception_to_r( ex );
    } catch(...) {
   ::Rf_error( "c++ exception (unknown reason)" );
    }
    return R_NilValue; // -Wall
}


(By the way, those **** in the signature are not really in the code - apparently the software censors submissions and decided to blank out the identifier type, which is an abbreviation of the phrase "symbolic expression" formed from the first letter of the first word and the first three letters of the second word.)

For full-rank cases this works as desired (and, by the way, I am quite impressed with the quality and flexibility of the Eigen library). I am not entirely sure that I have the reordering of the standard errors correct but that is a minor point which I will get to in the testing.

In the rank-deficient case, an assertion fails in the calculation of Rinv, reporting

Code: Select all

n = 22, p = 12, r = 11, df = 11
   -4.69042    -1.70561     -1.2792     -1.2792    -1.70561   -0.852803   -0.426401   -0.426401   -0.426401   -0.426401   -0.426401
   0.175734      2.2563  -0.0805823  -0.0805823    -1.28932    0.241747    0.564076   -0.322329    0.564076   -0.322329    0.564076
   0.175734   -0.273966     2.08738   -0.787044   -0.136877   -0.513289   -0.239535    0.684386   -0.239535   -0.273754    0.718605
   0.175734   -0.273966    0.112392     1.93331   -0.203507    -0.76315    0.678356  -0.0169589   -0.356137    0.627479   0.0339178
   0.175734    0.117265   0.0883996    0.112613    -1.83533   -0.344124  -0.0573539   -0.516185   0.0573539   -0.516185  -0.0573539
   0.175734    0.117265   0.0883996    0.112613    0.232613         1.5  -0.0833333  -0.0833333        0.75  -0.0833333  -0.0833333
   0.175734   -0.273966    0.112392    0.143177  -0.0267955   -0.378323    0.986013    0.140859   -0.253546   -0.535264   -0.366234
   0.175734   -0.273966    0.112392    0.143177  -0.0267955   -0.378323   0.0525165     -0.9759    -0.29277     0.29277     0.48795
   0.175734    0.117265   0.0883996    0.112613    0.232613   -0.454818   -0.130774   -0.204035    0.774597    0.258199   -0.258199
   0.175734    0.117265   0.0883996    0.112613    0.232613   -0.454818   -0.130774   -0.204035    0.111783    0.730297   -0.182574
   0.175734    0.117265   -0.336827   0.0141011   -0.159558 -0.00100954  -0.0342669   -0.261778   -0.221001   0.0945382    0.707107
   -0.213201    -0.161165    -0.136877    -0.203507     0.344124    -0.166667     0.112687     -0.09759 -1.06515e-16     0.273861     0.353553
           0     0.443203    0.0171096    0.0254383    -0.315447       -0.125    -0.295804 -7.53457e-17    -0.258199    -0.182574    -0.707107
           0            0      0.47907     0.195027   -0.0573539         0.25  2.65937e-17     0.341565     0.129099    -0.182574    -0.707107
           0            0            0     0.517246   -0.0573539         0.25    -0.338062    -0.048795    -0.129099     -0.63901    -0.353553
          -0           -0           -0           -0    -0.544862       -0.125   -0.0422577      0.29277     0.258199     -0.63901    -0.353553
           0            0            0            0            0     0.666667    0.0563436    -0.048795    -0.645497     0.365148   6.4763e-17
           0            0            0            0            0            0      1.01419     0.146385     0.387298     0.547723     0.707107
          -0           -0           -0           -0           -0           -0           -0      -1.0247    -0.387298     0.547723     0.707107
           0            0            0            0            0            0            0            0      1.29099    -0.456435     0.353553
           0            0            0            0            0            0            0            0            0      1.36931     0.353553
           0            0            0            0            0            0            0            0            0            0      1.41421
R: ../inst/include/Eigen/src/Core/Block.h:278: Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel, true>::Block(XprType&, Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel, true>::Index) [with XprType = Eigen::Matrix<double, -0x00000000000000001, 1>, int BlockRows = 1, int BlockCols = 1, bool InnerPanel = false, Eigen::Block<XprType, BlockRows, BlockCols, InnerPanel, true>::Index = long int]: Assertion `(i>=0) && ( ((BlockRows==1) && (BlockCols==XprType::ColsAtCompileTime) && i<xpr.rows()) ||((BlockRows==XprType::RowsAtCompileTime) && (BlockCols==1) && i<xpr.cols()))' failed.


I can provide the matrix A and the vector b that generates this example if that would help in creating a test case.

Am I doing the calculation of the standard errors, which are the diagonal elements of s^2 * Rinv * trans(Rinv), effectively? The matrix Rinv should be the inverse of the non-singular part of R (i.e. the first r rows and columns).

If this is a reasonable way of calculating these values, how can I avoid the assertion failure?

To generate Rinv I apply a solve method to an identity matrix. I looked around to see if there is an explicit triangular inverse (like Lapack's dtrtri) but didn't see one.

Thanks for creating Eigen and thanks for any advice.

Sorry, I misstated the definition of the standard errors. They are the square roots of the diagonal elements of s^2 * Rinv * trans(Rinv).

hi, the assertion is triggered by this line:

VectorXd se = mqr.colsPermutation() * Rinv.rowwise().norm() * s;

because mqr.colsPermutation() is a p x p permutation matrix while Rinv.rowwise().norm() is a rx1 column vector.

The pb only occurs when A is not full rank (r<p), but then the residual and standard error are zero anyway, so no need to compute them!

ggael wrote:hi, the assertion is triggered by this line:

VectorXd se = mqr.colsPermutation() * Rinv.rowwise().norm() * s;

because mqr.colsPermutation() is a p x p permutation matrix while Rinv.rowwise().norm() is a rx1 column vector.

Thank you. Dumb of me not to realize that.

ggael wrote:The pb only occurs when A is not full rank (r<p), but then the residual and standard error are zero anyway, so no need to compute them!

Not quite. One can still calculate the standard errors of the coefficients corresponding to the first r columns of A * mqr.colsPermutation(). Coefficients and standard errors for the remaining columns are not defined.

dmbates wrote:

ggael wrote:The pb only occurs when A is not full rank (r<p), but then the residual and standard error are zero anyway, so no need to compute them!

Not quite. One can still calculate the standard errors of the coefficients corresponding to the first r columns of A * mqr.colsPermutation(). Coefficients and standard errors for the remaining columns are not defined.

Of course, first I will need to look closely at the solve method to ensure that coefficients are solved correctly.

Hi GGael
Is there any inbuilt function that one can use to compute the inverse of a triangular matrix? There is the ubiquitous dtrtri from LAPACK which I can use but I would prefer to use something from Eigen if available.

Thanks

First of all make sure you really need to explicitly compute the inverse, it is often preferable to call solve multiple times. If so, then there is no shortcut, but you can still do something like:

Code: Select all

inv = mat.triangularView<Lower>().solve(MatrixXd::Identity(n,n));


Indeed, an inverse is the last resort that I would like to use. However, given that inversion can be performed recursively for a large matrix, does it make sense to use some form of recursion to compute the inverse?