如何找到矩形矩阵的列空间的基础？

Question

给定一个尺寸为 m x n（其中 n>m）的 numpy ndarray，我如何找到线性独立的列？

Answer 1

一种方法是使用LU decomposition。因子U 将与您的矩阵大小相同，但将是上三角形。在U 的每一行中，选择第一个非零元素：这些是枢轴元素，属于线性独立的列。一个独立的例子：

import numpy as np
from scipy.linalg import lu
A = np.array([[1, 2, 3], [2, 4, 2]])     # example for testing 
U = lu(A)[2]
lin_indep_columns = [np.flatnonzero(U[i, :])[0] for i in range(U.shape[0])]

输出：[0, 2]，表示 A 的第 0 列和第 2 列构成其列空间的基础。

Answer 2

@user6655984 的回答启发了这段代码，我在其中开发了一个函数而不是作者的最后一行代码（查找 U 的枢轴列），以便它可以处理更多样化的 A。

这里是：

import numpy as np
from scipy import linalg as LA

np.set_printoptions(precision=1, suppress=True)

A = np.array([[1, 4, 1, -1],
              [2, 5, 1, -2],
              [3, 6, 1, -3]])

P, L, U = LA.lu(A)

print('P', P, '', 'L', L, '', 'U', U, sep='\n')

输出：

P
[[0. 1. 0.]
 [0. 0. 1.]
 [1. 0. 0.]]

L
[[1.  0.  0. ]
 [0.3 1.  0. ]
 [0.7 0.5 1. ]]

U
[[ 3.   6.   1.  -3. ]
 [ 0.   2.   0.7 -0. ]
 [ 0.   0.  -0.  -0. ]]

我想出了这个功能：

def get_indices_for_linearly_independent_columns_of_A(U: np.ndarray) -> list:
    
    # I should first convert all "-0."s to "0." so that nonzero() can find them.
    U_copy = U.copy()
    U_copy[abs(U_copy) < 1.e-7] = 0

    # Because some rows in U may not have even one nonzero element,
    # I have to find the index for the first one in two steps.
    index_of_all_nonzero_cols_in_each_row = (
        [U_copy[i, :].nonzero()[0] for i in range(U_copy.shape[0])]
    )
    index_of_first_nonzero_col_in_each_row = (
        [indices[0] for indices in index_of_all_nonzero_cols_in_each_row
         if len(indices) > 0]
    )

    # Because two rows or more may have the same indices
    # for their first nonzero element, I should remove duplicates.
    unique_indices = sorted(list(set(index_of_first_nonzero_col_in_each_row)))
    return unique_indices

最后：

col_sp_A = A[:, get_indices_for_linearly_independent_columns_of_A(U)]
print(col_sp_A)

输出：

[[1 4]
 [2 5]
 [3 6]]

Answer 3

试试这个

def LU_decomposition(A):
    """
    Perform LU decompostion of a given matrix
    Args:
        A: the given matrix

    Returns: P, L and U, s.t. PA = LU

    """
    assert A.shape[0] == A.shape[1]
    N = A.shape[0]
    P_idx = np.arange(0, N, dtype=np.int16).reshape(-1, 1)
    for i in range(N - 1):
        pivot_loc = np.argmax(np.abs(A[i:, [i]])) + i
        if pivot_loc != i:
            A[[i, pivot_loc], :] = A[[pivot_loc, i], :]
            P_idx[[i, pivot_loc], :] = P_idx[[pivot_loc, i], :]
        A[i + 1:, i] /= A[i, i]
        A[i + 1:, i + 1:] -= A[i + 1:, [i]] * A[[i], i + 1:]
    U, L, P = np.zeros_like(A), np.identity(N), np.zeros((N, N), dtype=np.int16)
    for i in range(N):
        L[i, :i] = A[i, :i]
        U[i, i:] = A[i, i:]
        P[i, P_idx[i][0]] = 1
    return P.astype(np.float64), L, U



def get_bases(A):
    assert A.ndim == 2
    Q = gaussian_elimination(A)
    M, N = Q.shape
    pivot_idxs = []
    for i in range(M):
        j = i
        while j < N and abs(Q[i, j]) < 1e-5:
            j += 1
        if j < N:
            pivot_idxs.append(j)
    return A[:, list(set(pivot_idxs))]