等宽矩阵(a)相乘a %*% x = b的逆运算solve(a,b)=x

> a <- matrix(1:12,2,6)
> x <- matrix(1:12,6,2)
> a
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    3    5    7    9   11
[2,]    2    4    6    8   10   12
> x
     [,1] [,2]
[1,]    1    7
[2,]    2    8
[3,]    3    9
[4,]    4   10
[5,]    5   11
[6,]    6   12
> b <- a %*% x
> b
     [,1] [,2]
[1,]  161  377
[2,]  182  434
> solve(a,b)
Error in solve.default(a, b) : 'a' (2 x 6) must be square

   F77_CALL(dgesv)(&n, &p, avals, &n, ipiv, REAL(B), &n, &info);
    if (info < 0)
        error(_("argument %d of Lapack routine %s had invalid value"),
              -info, "dgesv");
    if (info > 0)
        error(_("Lapack routine %s: system is exactly singular: U[%d,%d] = 0"),
              "dgesv", info, info);

> a <- matrix(1:16,4,4)
> solve(a, a %*% a)
Error in solve.default(a, a %*% a) : 
  Lapack routine dgesv: system is exactly singular: U[3,3] = 0

> a <- matrix(1:4,2,2)
> a
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> solve(a)   # solve(a) , 不给b的话, 其实b 默认就是和a维度一致并且对角线为1的矩阵.
     [,1] [,2]
[1,]   -2  1.5
[2,]    1 -0.5
> solve(a, diag(1,2,2))   # 可以看到
     [,1] [,2]
[1,]   -2  1.5
[2,]    1 -0.5
> diag(1,2,2)
     [,1] [,2]
[1,]    1    0
[2,]    0    1
> a %*% solve(a)   #  因为solve(a,b) = x, a %*% x = b, 所以 a %*% solve(a) = b = diag(1,2,2)
     [,1] [,2]
[1,]    1    0
[2,]    0    1

> solve(a, c(2,3))
[1] 0.5 0.5
> a %*% solve(a, c(2,3))
     [,1]
[1,]    2
[2,]    3

solve                   package:base                   R Documentation

Solve a System of Equations

Description:

     This generic function solves the equation ‘a %*% x = b’ for ‘x’,
     where ‘b’ can be either a vector or a matrix.

Usage:

     solve(a, b, ...)
     
     ## Default S3 method:
     solve(a, b, tol, LINPACK = FALSE, ...)
     
Arguments:

       a: a square numeric or complex matrix containing the
          coefficients of the linear system.  Logical matrices are
          coerced to numeric.

       b: a numeric or complex vector or matrix giving the right-hand
          side(s) of the linear system.  If missing, ‘b’ is taken to be
          an identity matrix and ‘solve’ will return the inverse of
          ‘a’.

     tol: the tolerance for detecting linear dependencies in the
          columns of ‘a’.  The default is ‘.Machine$double.eps’. Not
          currently used with complex matrices ‘a’.

 LINPACK: logical.  Defunct and ignored.

     ...: further arguments passed to or from other methods

/* Real case of solve.default */
static SEXP La_solve(SEXP A, SEXP Bin, SEXP tolin)
{
    int n, p;
    double *avals, anorm, rcond, tol = asReal(tolin), *work;
    SEXP B, Adn, Bdn;

    if (!(isMatrix(A) && 
          (TYPEOF(A) == REALSXP || TYPEOF(A) == INTSXP || TYPEOF(A) == LGLSXP)))
        error(_("'a' must be a numeric matrix"));
    int *Adims = INTEGER(coerceVector(getAttrib(A, R_DimSymbol), INTSXP));
    n = Adims[0];
    if(n == 0) error(_("'a' is 0-diml"));
    int n2 = Adims[1];
    if(n2 != n) error(_("'a' (%d x %d) must be square"), n, n2);
    Adn = getAttrib(A, R_DimNamesSymbol);

    if (isMatrix(Bin)) {
        int *Bdims = INTEGER(coerceVector(getAttrib(Bin, R_DimSymbol), INTSXP));
        p = Bdims[1];
        if(p == 0) error(_("no right-hand side in 'b'"));
        int p2 = Bdims[0];
        if(p2 != n)
            error(_("'b' (%d x %d) must be compatible with 'a' (%d x %d)"),
                  p2, p, n, n);
        PROTECT(B = allocMatrix(REALSXP, n, p));
        SEXP Bindn =  getAttrib(Bin, R_DimNamesSymbol);
        // This is somewhat odd, but Matrix relies on dropping NULL dimnames
        if (!isNull(Adn) || !isNull(Bindn)) {
            // rownames(ans) = colnames(A), colnames(ans) = colnames(Bin)
            Bdn = allocVector(VECSXP, 2);
            if (!isNull(Adn)) SET_VECTOR_ELT(Bdn, 0, VECTOR_ELT(Adn, 1));
            if (!isNull(Bindn)) SET_VECTOR_ELT(Bdn, 1, VECTOR_ELT(Bindn, 1));
            if (!isNull(VECTOR_ELT(Bdn, 0)) || !isNull(VECTOR_ELT(Bdn, 1)))
                setAttrib(B, R_DimNamesSymbol, Bdn);
        }
    } else {
        p = 1;
        if(length(Bin) != n)
            error(_("'b' (%d x %d) must be compatible with 'a' (%d x %d)"),
                  length(Bin), p, n, n);        
        PROTECT(B = allocVector(REALSXP, n));
        if (!isNull(Adn)) setAttrib(B, R_NamesSymbol, VECTOR_ELT(Adn, 1));
    }
    PROTECT(Bin = coerceVector(Bin, REALSXP));
    Memcpy(REAL(B), REAL(Bin), (size_t)n * p);
    
    int *ipiv = (int *) R_alloc(n, sizeof(int));

    /* work on a copy of A */
    if (!isReal(A)) {
        A = coerceVector(A, REALSXP);
        avals = REAL(A);
    } else {
        avals = (double *) R_alloc((size_t)n * n, sizeof(double));
        Memcpy(avals, REAL(A), (size_t)n * n);
    }
    PROTECT(A);
    int info;
    F77_CALL(dgesv)(&n, &p, avals, &n, ipiv, REAL(B), &n, &info);
    if (info < 0)
        error(_("argument %d of Lapack routine %s had invalid value"),
              -info, "dgesv");
    if (info > 0)
        error(_("Lapack routine %s: system is exactly singular: U[%d,%d] = 0"),
              "dgesv", info, info);
    if(tol > 0) {
        char one[2] = "1";
        anorm = F77_CALL(dlange)(one, &n, &n, REAL(A), &n, (double*) NULL);
        work = (double *) R_alloc(4*(size_t)n, sizeof(double));
        F77_CALL(dgecon)(one, &n, avals, &n, &anorm, &rcond, work, ipiv, &info);
        if (rcond < tol)
            error(_("system is computationally singular: reciprocal condition number = %g"),
                  rcond);
    }
    UNPROTECT(3); /* B, Bin, A */
    return B;
}