CODEWITHSUNDEEP
numpy
Demetri Pananos
Demetri Pananos

Reputation: 7404

How can I find a basis for the column space of a rectangular matrix?

Given a numpy ndarray with dimensions m by n (where n>m), how can I find the linearly independent columns?

Upvotes: 1

Views: 3466

Answers (3)

Rancho Xia
Rancho Xia

Reputation: 1

Try this one



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))]

Upvotes: 0

Shahrokh Bah
Shahrokh Bah

Reputation: 410

@user6655984's answer inspired this code, where I developed a function instead of the author's last line of code (finding pivot columns of U) so that it can handle more diverse A's.

Here it is:

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')

Output:

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. ]]

I came up with this function:

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

Finally:

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

Output:

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

Upvotes: 0

user6655984
user6655984

Reputation:

One way is to use the LU decomposition. The factor U will be of the same size as your matrix, but will be upper-triangular. In each row of U, pick the first nonzero element: these are pivot elements, which belong to linearly independent columns. A self-contained example: 

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])]

Output: [0, 2], which means the 0th and 2nd columns of A form a basis for its column space.

Upvotes: 3

Related Questions