rpca.py

import numpy as np

def shrink(X, tau):
    return np.sign(X) * np.maximum(np.abs(X) - tau, 0)

def ialm(D, threshold=1e-3):
    """
    Performs low-rank and sparse decomposition of the input matrix D using
    the Inexact Augmented Lagrange Multiplier (IALM) method.
    
    The algorithm decomposes D into:
        D = A + E
    where A is a low-rank matrix and E is a sparse error matrix.
    
    Parameters:
        D (np.ndarray): The input data matrix of shape (m, n).
        threshold (float): Convergence threshold based on the residual change.
        
    Returns:
        A (np.ndarray): The recovered low-rank component.
        E (np.ndarray): The recovered sparse error component.
        iterations (int): The number of iterations executed.
    """
    m, n = D.shape
    iterations = 0
    Y = np.zeros((m, n))
    A = np.zeros((m, n))
    E = np.zeros((m, n))

    # TODO: Fill this functions
    # Important Notes: You can use any functions in numpy to implement svd!
    # Recommend you to use frobenius norm in numpy

    mu = 1
    p = 1.1
    mu_max = mu * 1e-7
    lam = 1

    norm_D = np.linalg.norm(D, 'fro')

    for i in range(1000):
        iterations += 1

        Temp_A = D - E + (1.0 / mu) * Y
        U, S, Vt = np.linalg.svd(Temp_A, full_matrices=False)
        s_thresholded = np.diag(shrink(S, 1.0/mu))
        A_next = (U @ s_thresholded) @ Vt

        Temp_E = D - A_next + (1.0 / mu) * Y
        E_next = shrink(Temp_E, lam / mu)

        R = D - A_next - E_next  # residual
        Y_next = Y + mu * R

        err = np.linalg.norm(R, 'fro') / norm_D
        if err < threshold:
            # A, E, Y = A_next, E_next, Y_next
            break

        A, E, Y = A_next, E_next, Y_next

        mu = min(p * mu, mu_max)

    print(iterations)
    return A, E, iterations