Skip to content

project_out_subspace

project_out_subspace(A, B, M=None)

Performs a least squares projection on subspace A so that it does not span space B. Finds C as close as possible to A such that C is orthogonal to B.

    argmin_C ||C - A||^2_F st. C^T B = 0
    ||C - A||^2_F st. C^T B = 0
    C'MC - 2 C'MA + mu' C'B

    [ M  B ] [ C ]  = [ MA ]
    [ B' 0 ] [ mu ]   [ 0 ]

    C = M^-1 (MA - B mu)
    B' C = 0 -> B' M^-1 (MA - B mu) = 0 ->  mu = ( B' M^-1 B) B' A
    C = A - M^-1 B  (B' M^-1 B)^(-1) B' A

Parameters:

Name Type Description Default
A (n, m) float numpy array

Subspace of interes.

 required
B (n, k) float numpy array

Subspace to project out. We want C to be orthogonal to B

 required
M (n, n) float numpy array

Mass matrix defining the metric for projection. If None, set to identity matrix

 None

Returns:

Name Type Description
C (n, m) float numpy array

Projection of A
Source code in src\fast_cody\project_out_subspace.py 
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
def project_out_subspace(A, B, M=None):
    """ Performs a least squares projection on subspace A so that it does not span space B.
    Finds C as close as possible to A such that C is orthogonal to B.
    ```
        argmin_C ||C - A||^2_F st. C^T B = 0
        ||C - A||^2_F st. C^T B = 0
        C'MC - 2 C'MA + mu' C'B

        [ M  B ] [ C ]  = [ MA ]
        [ B' 0 ] [ mu ]   [ 0 ]

        C = M^-1 (MA - B mu)
        B' C = 0 -> B' M^-1 (MA - B mu) = 0 ->  mu = ( B' M^-1 B) B' A
        C = A - M^-1 B  (B' M^-1 B)^(-1) B' A
    ```

    Parameters
    ----------
    A : (n, m) float numpy array
        Subspace of interes.
    B : (n, k) float numpy array
        Subspace to project out. We want C to be orthogonal to B
    M : (n, n) float numpy array
        Mass matrix defining the metric for projection. If None, set to identity matrix

    Returns
    -------
    C : (n, m) float numpy array
        Projection of A

    """


    assert(B.shape[0] == A.shape[0])
    if M is None:
        M = sp.sparse.identity(A.shape[0]).tocsc()

    Z = sp.sparse.csc_matrix((B.shape[1], B.shape[1]))
    Bsp = sp.sparse.csc_matrix(B)
    Q = vstack((hstack([M, Bsp]), hstack((Bsp.T, Z))))
    # Q = vstack(hstack((M, B)), hstack((B.T, Z)))

    z = sp.sparse.csc_matrix((B.shape[1], A.shape[1]))

    Asp = sp.sparse.csc_matrix(A)
    rhs = vstack(( M @ Asp, z))

    Cmu = umfpack_lu_solve(Q, rhs.todense())
    #sp.sparse.linalg.spsolve(Q, rhs)
    C = Cmu[:A.shape[0], :]
    # Cd = C.toarray()
    return C