[R] off topic: C/C++ codes for pseudo inverse

[Torsten Hothorn]

On Wed, 16 Jun 2004, Mahbub Latif wrote:

> Hi,
>
> I am looking for C/C++ codes for computing generalized
> inverse of a matrix. Can anyone help me in this
> regard?
>

experimental, undocumented and for square matrices only:

SEXP svd (SEXP x) {

    /* experimental */

    SEXP jobu, jobv, u, v, method, dummy, ans;
    int i, p;

    if (!isMatrix(x) || !isReal(x))
        error("x is not a real matrix");

    PROTECT(method = allocVector(STRSXP, 1));
    SET_STRING_ELT(method, 0, mkChar("dgesdd"));
    PROTECT(jobu = allocVector(STRSXP, 1));
    SET_STRING_ELT(jobu, 0, mkChar("S"));
    PROTECT(jobv = allocVector(STRSXP, 1));
    SET_STRING_ELT(jobv, 0, mkChar(""));
    p = INTEGER(getAttrib(x, R_DimSymbol))[0];
    PROTECT(u = allocMatrix(REALSXP, p, p));
    PROTECT(v = allocMatrix(REALSXP, p, p));
    for (i = 0; i < p*p; i++) {
        REAL(u)[i] = 0;
        REAL(v)[i] = 0;
    }
    PROTECT(dummy = allocVector(REALSXP, p));
    PROTECT(ans = La_svd(jobu, jobv, x, dummy, u, v, method));
    UNPROTECT(7);
    return(ans);
}

void cMPinv (SEXP x, double tol, double *ans, double *rank) {

    SEXP svdx, d, u, vt;
    int i, j, p, k, *positive;
    double *dd, *du, *dvt;

    PROTECT(svdx = svd(x));
    d = VECTOR_ELT(svdx, 0);
    dd = REAL(d);
    u = VECTOR_ELT(svdx, 1);
    du = REAL(u);
    vt = VECTOR_ELT(svdx, 2);
    dvt = REAL(vt);
    p = LENGTH(d);

    tol = tol * dd[0];
    if (tol < 0.0) tol = 0.0;

    positive = Calloc(p, int);

    rank[0] = 0.0;
    for (i = 0; i < p; i++) {
        if (dd[i] > tol) {
            positive[i] = 1;
            rank[0]++;
        }
    }

    for (j = 0; j < p; j++) {
        if (positive[j]) {
            for (i = 0; i < p; i++)
                du[j * p + i] *= (1 / dd[j]);
        }
    }

    for (i = 0; i < p; i++) {
        for (j = 0; j < p; j++) {
            ans[j * p + i] = 0.0;
            for (k = 0; k < p; k++) {
                if (positive[k])
                    ans[j * p + i] += dvt[i * p + k] * du[p * k + j];
            }
        }
    }

    Free(positive);
    UNPROTECT(1);
}

SEXP MPinv (SEXP x, SEXP tol) {

    SEXP ans, mp, rank;
    int p;

    if (!isMatrix(x) || !isReal(x))
        error("x is not a real matrix");

    if (INTEGER(getAttrib(x, R_DimSymbol))[0] !=
        INTEGER(getAttrib(x, R_DimSymbol))[1])
        error("x is not a square matrix");

    if (!isReal(tol) || LENGTH(tol) != 1)
        error("tol is not a scalar real");

    p = INTEGER(getAttrib(x, R_DimSymbol))[0];

    PROTECT(ans = allocVector(VECSXP, 2));
    SET_VECTOR_ELT(ans, 0, mp = allocMatrix(REALSXP, p, p));
    SET_VECTOR_ELT(ans, 1, rank = allocVector(REALSXP, 1));

    cMPinv(x, REAL(tol)[0], REAL(mp), REAL(rank));

    UNPROTECT(1);
    return(ans);
}

you need to include `La_svd' from R.

Torsten

> Thanks,
>
> Mahbub.
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>
>