为 R 中的任意 m * n 矩阵 A 求解齐次系统 Ax = 0（找到 A 的空 space 基）
Solve homogenous system Ax = 0 for any m * n matrix A in R (find null space basis for A)
当 A 是 R 中的任何 m * n 矩阵（不一定是平方矩阵）时，如何求解齐次系统 Ax = 0？
# A=[-0.1 0.1]= 1x2 matrix; x=2x1 to be found; 0: 1x1 zero matrix
A <- t(matrix(c(-0.1,0.1)))
这道题好像等价于找到一个Rn -> Rm（不能做上标；sorry）线性变换的核（null space）。
Anyway, the solution for the above specific matrix A will suffice to me.
我们可以eye-spot它，x = (a, a)，其中a是任意值。
经典/text-book解决方案
下面的函数 NullSpace 使用上述理论找到 A 的 null space。在情况 1 中，null space 平凡地为零；而在情况 2 到 4 中返回一个矩阵，其列跨越 null space.
NullSpace <- function (A) {
m <- dim(A)[1]; n <- dim(A)[2]
## QR factorization and rank detection
QR <- base::qr.default(A)
r <- QR$rank
## cases 2 to 4
if ((r < min(m, n)) || (m < n)) {
R <- QR$qr[1:r, , drop = FALSE]
P <- QR$pivot
F <- R[, (r + 1):n, drop = FALSE]
I <- base::diag(1, n - r)
B <- -1.0 * base::backsolve(R, F, r)
Y <- base::rbind(B, I)
X <- Y[base::order(P), , drop = FALSE]
return(X)
}
## case 1
return(base::matrix(0, n, 1))
}
使用您的示例矩阵，它正确 returns 空值 space。
A1 <- matrix(c(-0.1, 0.1), 1, 2)
NullSpace(A1)
# [,1]
#[1,] 1
#[2,] 1
我们也可以尝试一个随机的例子。
set.seed(0)
A2 <- matrix(runif(10), 2, 5)
# [,1] [,2] [,3] [,4] [,5]
#[1,] 0.8966972 0.3721239 0.9082078 0.8983897 0.6607978
#[2,] 0.2655087 0.5728534 0.2016819 0.9446753 0.6291140
X <- NullSpace(A2)
# [,1] [,2] [,3]
#[1,] -1.0731435 -0.393154 -0.3481344
#[2,] 0.1453199 -1.466849 -0.9368564
#[3,] 1.0000000 0.000000 0.0000000
#[4,] 0.0000000 1.000000 0.0000000
#[5,] 0.0000000 0.000000 1.0000000
## X indeed solves A2 %*% X = 0
A2 %*% X
# [,1] [,2] [,3]
#[1,] 2.220446e-16 -1.110223e-16 -2.220446e-16
#[2,] 5.551115e-17 -1.110223e-16 -1.110223e-16
寻找正交基
我的函数 NullSpace returns 一个 non-orthogonal 基础为 null space。另一种方法是将 QR 因式分解应用于 t(A)（A 的转置），获得等级 r，并保留最后的 (n - r) 列Q矩阵。这给出了空值 space.
的 正交基
pracma 包中的 nullspace 函数是一个现有的实现。我们拿上面的矩阵A2来演示一下。
library(pracma)
X2 <- nullspace(A2)
# [,1] [,2] [,3]
#[1,] -0.67453687 -0.24622524 -0.2182437
#[2,] 0.27206765 -0.69479881 -0.4260258
#[3,] 0.67857488 0.07429112 0.0200459
#[4,] -0.07098962 0.62990141 -0.2457700
#[5,] -0.07399872 -0.23309397 0.8426547
## it indeed solves A2 %*% X = 0
A2 %*% X2
# [,1] [,2] [,3]
#[1,] 2.567391e-16 1.942890e-16 0.000000e+00
#[2,] 6.938894e-17 -5.551115e-17 -1.110223e-16
## it has orthonormal columns
round(crossprod(X2), 15)
# [,1] [,2] [,3]
#[1,] 1 0 0
#[2,] 0 1 0
#[3,] 0 0 1
附录：图片的 Markdown（需要 MathJax 支持）
Let $A$ be $m \times n$, then its null space is $\{x: Ax = 0\}$. To find a solution of $Ax = 0$, the conventional method is Gaussian elimination that reduces $A$ into a row echelon form. However, let's consider the QR factorization (with column pivoting) approach, where a sequence of Householder transforms are applied to both sides of $Ax = 0$, reducing the equation to $RP'x = 0$, where $P$ is an $n \times n$ column permutation matrix. What $R$ looks like depends on the relationship of $m$ and $n$, as well as the rank of $A$, denoted by $r$.
1. If $m \ge n = r$, $R$ is an $n \times n$ full-rank upper triangular matrix, which looks like $$\begin{pmatrix} \times & \times & \times & \times \ & \times & \times & \times \ & & \times & \times \ & & & \times\end{pmatrix}$$
2. If $m \ge n > r$, $R$ is an $n \times n$ rank-deficient upper triangular matrix, which looks like $$\begin{pmatrix} \times & \times & \times & \times \ & \times & \times & \times \ & & \times & \times \ & & & 0\end{pmatrix}$$
3. If $r = m < n$, $R$ is an $m \times n$ full-rank matrix which looks like $$\begin{pmatrix} \times & \times & \times & \times & \times & \times & \times \ & \times & \times & \times & \times & \times & \times\ & & \times & \times & \times & \times & \times\ & & & \times & \times & \times & \times\end{pmatrix}$$
4. If $r < m < n$, $R$ is an $m \times n$ rank-deficient matrix which looks like $$\begin{pmatrix} \times & \times & \times & \times & \times & \times & \times \ & \times & \times & \times & \times & \times & \times\ & & \times & \times & \times & \times & \times\ & & & 0 & 0 & 0 & 0\end{pmatrix}$$.
In all cases, the first $r$ non-zero rows of $R$ can be partiontioned into $\begin{pmatrix} U & F\end{pmatrix}$, where $U$ is an $r \times r$ full-rank upper triangular matrix and $F$ is an $r \times (n - r)$ rectangular matrix. The null space of $A$ has dimension $(n - r)$ and can be characterized by an $ n \times (n - r)$ matrix $X$, such that $RP'X = 0$. In practice, $X$ is obtained in two steps.
1. Let $Y$ be an $ n \times (n - r)$ matrix and solve $RY = 0$. Clearly $Y$ can not be uniquely determined as the linear system has $(n - r)$ free parameters. A common solution is to find $Y = \left(\begin{smallmatrix} B \ I \end{smallmatrix}\right)$, where $I$ is an $(n - r) \times (n - r)$ identity matrix. Then $B$ can be uniquely solved from $UB = -F$ using back substitution.
2. Solve $P'X = Y$ for $X = PY$, which is simply a row permutation of $Y$.
当 A 是 R 中的任何 m * n 矩阵（不一定是平方矩阵）时，如何求解齐次系统 Ax = 0？
# A=[-0.1 0.1]= 1x2 matrix; x=2x1 to be found; 0: 1x1 zero matrix
A <- t(matrix(c(-0.1,0.1)))
这道题好像等价于找到一个Rn -> Rm（不能做上标；sorry）线性变换的核（null space）。
Anyway, the solution for the above specific matrix
Awill suffice to me.
我们可以eye-spot它，x = (a, a)，其中a是任意值。
经典/text-book解决方案
下面的函数 NullSpace 使用上述理论找到 A 的 null space。在情况 1 中，null space 平凡地为零；而在情况 2 到 4 中返回一个矩阵，其列跨越 null space.
NullSpace <- function (A) {
m <- dim(A)[1]; n <- dim(A)[2]
## QR factorization and rank detection
QR <- base::qr.default(A)
r <- QR$rank
## cases 2 to 4
if ((r < min(m, n)) || (m < n)) {
R <- QR$qr[1:r, , drop = FALSE]
P <- QR$pivot
F <- R[, (r + 1):n, drop = FALSE]
I <- base::diag(1, n - r)
B <- -1.0 * base::backsolve(R, F, r)
Y <- base::rbind(B, I)
X <- Y[base::order(P), , drop = FALSE]
return(X)
}
## case 1
return(base::matrix(0, n, 1))
}
使用您的示例矩阵，它正确 returns 空值 space。
A1 <- matrix(c(-0.1, 0.1), 1, 2)
NullSpace(A1)
# [,1]
#[1,] 1
#[2,] 1
我们也可以尝试一个随机的例子。
set.seed(0)
A2 <- matrix(runif(10), 2, 5)
# [,1] [,2] [,3] [,4] [,5]
#[1,] 0.8966972 0.3721239 0.9082078 0.8983897 0.6607978
#[2,] 0.2655087 0.5728534 0.2016819 0.9446753 0.6291140
X <- NullSpace(A2)
# [,1] [,2] [,3]
#[1,] -1.0731435 -0.393154 -0.3481344
#[2,] 0.1453199 -1.466849 -0.9368564
#[3,] 1.0000000 0.000000 0.0000000
#[4,] 0.0000000 1.000000 0.0000000
#[5,] 0.0000000 0.000000 1.0000000
## X indeed solves A2 %*% X = 0
A2 %*% X
# [,1] [,2] [,3]
#[1,] 2.220446e-16 -1.110223e-16 -2.220446e-16
#[2,] 5.551115e-17 -1.110223e-16 -1.110223e-16
寻找正交基
我的函数 NullSpace returns 一个 non-orthogonal 基础为 null space。另一种方法是将 QR 因式分解应用于 t(A)（A 的转置），获得等级 r，并保留最后的 (n - r) 列Q矩阵。这给出了空值 space.
pracma 包中的 nullspace 函数是一个现有的实现。我们拿上面的矩阵A2来演示一下。
library(pracma)
X2 <- nullspace(A2)
# [,1] [,2] [,3]
#[1,] -0.67453687 -0.24622524 -0.2182437
#[2,] 0.27206765 -0.69479881 -0.4260258
#[3,] 0.67857488 0.07429112 0.0200459
#[4,] -0.07098962 0.62990141 -0.2457700
#[5,] -0.07399872 -0.23309397 0.8426547
## it indeed solves A2 %*% X = 0
A2 %*% X2
# [,1] [,2] [,3]
#[1,] 2.567391e-16 1.942890e-16 0.000000e+00
#[2,] 6.938894e-17 -5.551115e-17 -1.110223e-16
## it has orthonormal columns
round(crossprod(X2), 15)
# [,1] [,2] [,3]
#[1,] 1 0 0
#[2,] 0 1 0
#[3,] 0 0 1
附录：图片的 Markdown（需要 MathJax 支持）
Let $A$ be $m \times n$, then its null space is $\{x: Ax = 0\}$. To find a solution of $Ax = 0$, the conventional method is Gaussian elimination that reduces $A$ into a row echelon form. However, let's consider the QR factorization (with column pivoting) approach, where a sequence of Householder transforms are applied to both sides of $Ax = 0$, reducing the equation to $RP'x = 0$, where $P$ is an $n \times n$ column permutation matrix. What $R$ looks like depends on the relationship of $m$ and $n$, as well as the rank of $A$, denoted by $r$.
1. If $m \ge n = r$, $R$ is an $n \times n$ full-rank upper triangular matrix, which looks like $$\begin{pmatrix} \times & \times & \times & \times \ & \times & \times & \times \ & & \times & \times \ & & & \times\end{pmatrix}$$
2. If $m \ge n > r$, $R$ is an $n \times n$ rank-deficient upper triangular matrix, which looks like $$\begin{pmatrix} \times & \times & \times & \times \ & \times & \times & \times \ & & \times & \times \ & & & 0\end{pmatrix}$$
3. If $r = m < n$, $R$ is an $m \times n$ full-rank matrix which looks like $$\begin{pmatrix} \times & \times & \times & \times & \times & \times & \times \ & \times & \times & \times & \times & \times & \times\ & & \times & \times & \times & \times & \times\ & & & \times & \times & \times & \times\end{pmatrix}$$
4. If $r < m < n$, $R$ is an $m \times n$ rank-deficient matrix which looks like $$\begin{pmatrix} \times & \times & \times & \times & \times & \times & \times \ & \times & \times & \times & \times & \times & \times\ & & \times & \times & \times & \times & \times\ & & & 0 & 0 & 0 & 0\end{pmatrix}$$.
In all cases, the first $r$ non-zero rows of $R$ can be partiontioned into $\begin{pmatrix} U & F\end{pmatrix}$, where $U$ is an $r \times r$ full-rank upper triangular matrix and $F$ is an $r \times (n - r)$ rectangular matrix. The null space of $A$ has dimension $(n - r)$ and can be characterized by an $ n \times (n - r)$ matrix $X$, such that $RP'X = 0$. In practice, $X$ is obtained in two steps.
1. Let $Y$ be an $ n \times (n - r)$ matrix and solve $RY = 0$. Clearly $Y$ can not be uniquely determined as the linear system has $(n - r)$ free parameters. A common solution is to find $Y = \left(\begin{smallmatrix} B \ I \end{smallmatrix}\right)$, where $I$ is an $(n - r) \times (n - r)$ identity matrix. Then $B$ can be uniquely solved from $UB = -F$ using back substitution.
2. Solve $P'X = Y$ for $X = PY$, which is simply a row permutation of $Y$.