Python scipy.linalg.orth() Examples

def random_walk(G,initial_prob,subspace_dim=3,walk_steps=3):
    """
    Start a random walk with probability distribution p_initial. 
    Transition matrix needs to be calculated according to adjacent matrix G.
    
    """
    assert type(initial_prob) == np.ndarray, "Initial probability distribution is \
                                             not a numpy array"
       
    # Transform the adjacent matrix to a laplacian matrix P
    P = adj_to_Laplacian(G)
    
    Prob_Matrix = np.zeros((G.shape[0], subspace_dim))
    Prob_Matrix[:,0] = initial_prob
    for i in range(1,subspace_dim):
        Prob_Matrix[:,i] = np.dot(Prob_Matrix[:,i-1], P)
     
    Orth_Prob_Matrix = splin.orth(Prob_Matrix)
    
    for i in range(walk_steps):
        temp = np.dot(Orth_Prob_Matrix.T, P)
        Orth_Prob_Matrix = splin.orth(temp.T)
    
    return Orth_Prob_Matrix 

def compare_solutions(A,B,m):
    n = A.shape[0]

    numpy.random.seed(0)

    V = rand(n,m)
    X = linalg.orth(V)

    eigs,vecs = lobpcg(A, X, B=B, tol=1e-5, maxiter=30)
    eigs.sort()

    #w,v = symeig(A,B)
    w,v = eig(A,b=B)
    w.sort()

    assert_almost_equal(w[:int(m/2)],eigs[:int(m/2)],decimal=2)

    #from pylab import plot, show, legend, xlabel, ylabel
    #plot(arange(0,len(w[:m])),w[:m],'bx',label='Results by symeig')
    #plot(arange(0,len(eigs)),eigs,'r+',label='Results by lobpcg')
    #legend()
    #xlabel(r'Eigenvalue $i$')
    #ylabel(r'$\lambda_i$')
    #show() 

def find_shared_subspace(mat1, mat2, sin_thres=0.05, cos_thres=1 / math.sqrt(2), mat2_vec=False,
                         assume_orthonomal=False, get_angle=True):
    if mat2_vec:
        mat2 = mat2[:, None]
    if not assume_orthonomal:
        mat1 = orth(mat1)
        mat2 = orth(mat2)
    cross_prod = np.dot(mat1.T, mat2)
    singular = np.linalg.svd(cross_prod)
    shared = sum(singular[1] > sin_thres)
    if not get_angle:
        return None, shared
    costheta = min(singular[1])
    if costheta < cos_thres:
        theta = math.acos(min(1, costheta))
    else:
        if mat1.shape[1] < mat2.shape[1]:
            sintheta = np.linalg.norm(x=mat1 - np.dot(mat2, cross_prod.T), ord=2)
        else:
            sintheta = np.linalg.norm(x=mat2.T - np.dot(mat1, cross_prod), ord=2)
        theta = math.asin(min(1, sintheta))
    return 180 * theta / math.pi, shared 

def __random_walk(G, initial_prob, subspace_dim=3, walk_steps=3):
    """
    Start a random walk with probability distribution p_initial.
    Transition matrix needs to be calculated according to adjacent matrix G.

    """
    assert type(initial_prob) == np.ndarray, "Initial probability distribution is \
                                             not a numpy array"

    # Transform the adjacent matrix to a laplacian matrix P
    P = __adj_to_Laplacian(G)

    Prob_Matrix = np.zeros((G.shape[0], subspace_dim))
    Prob_Matrix[:, 0] = initial_prob
    for i in range(1, subspace_dim):
        Prob_Matrix[:, i] = np.dot(Prob_Matrix[:, i - 1], P)

    Orth_Prob_Matrix = splin.orth(Prob_Matrix)

    for i in range(walk_steps):
        temp = np.dot(Orth_Prob_Matrix.T, P)
        Orth_Prob_Matrix = splin.orth(temp.T)

    return Orth_Prob_Matrix 

def test_fit_elbows():
    n=10; elbows=3
    np.random.seed(1)
    x = np.random.binomial(1, 0.6, (n ** 2)).reshape(n, n)
    xorth = orth(x)
    d = np.zeros(xorth.shape[0])
    for i in range(0, len(d), int(len(d) / (elbows + 1))):
        d[:i] += 10
    X = xorth.T.dot(np.diag(d)).dot(xorth)

    Xs = [X, X]

    ajive = AJIVE(n_elbows=2)
    ajive = ajive.fit(Xs)

    np.testing.assert_equal(list(ajive.init_signal_ranks_.values())[0], 4) 

def compare_solutions(A,B,m):
    n = A.shape[0]

    np.random.seed(0)

    V = rand(n,m)
    X = linalg.orth(V)

    eigs,vecs = lobpcg(A, X, B=B, tol=1e-5, maxiter=30)
    eigs.sort()

    w,v = eig(A,b=B)
    w.sort()

    assert_almost_equal(w[:int(m/2)],eigs[:int(m/2)],decimal=2) 

def generate_data(n=10, elbows=3, seed=1):
    """
    Generate data matrix with a specific number of elbows on scree plot
    """
    np.random.seed(seed)
    x = np.random.binomial(1, 0.6, (n ** 2)).reshape(n, n)
    xorth = orth(x)
    d = np.zeros(xorth.shape[0])
    for i in range(0, len(d), int(len(d) / (elbows + 1))):
        d[:i] += 10
    A = xorth.T.dot(np.diag(d)).dot(xorth)
    return A, d 

def generate_data(n=10, elbows=3, seed=1):
    """
    Generate data matrix with a specific number of elbows on scree plot
    """
    np.random.seed(seed)
    x = np.random.binomial(1, 0.6, (n ** 2)).reshape(n, n)
    xorth = orth(x)
    d = np.zeros(xorth.shape[0])
    for i in range(0, len(d), int(len(d) / (elbows + 1))):
        d[:i] += 10
    A = xorth.T.dot(np.diag(d)).dot(xorth)
    return A, d 

def generate_data(n=10, elbows=3, seed=1):
    """
    Generate data matrix with a specific number of elbows on scree plot
    """
    np.random.seed(seed)
    x = np.random.binomial(1, 0.6, (n ** 2)).reshape(n, n)
    xorth = orth(x)
    d = np.zeros(xorth.shape[0])
    for i in range(0, len(d), int(len(d) / (elbows + 1))):
        d[:i] += 10
    A = xorth.T.dot(np.diag(d)).dot(xorth)
    return A, d 

def gen_union_of_subspaces(ambient_dim, subspace_dim, num_subspaces, num_points_per_subspace, noise_level=0.0):
    """This funtion generates a union of subspaces under random model, i.e., 
    subspaces are independently and uniformly distributed in the ambient space,
    data points are independently and uniformly distributed on the unit sphere of each subspace

    Parameters
    -----------
    ambient_dim : int
        Dimention of the ambient space
    subspace_dim : int
        Dimension of each subspace (all subspaces have the same dimension)
    num_subspaces : int
        Number of subspaces to be generated
    num_points_per_subspace : int
        Number of data points from each of the subspaces
    noise_level : float
        Amount of Gaussian noise on data
		
    Returns
    -------
    data : shape (num_subspaces * num_points_per_subspace) by ambient_dim
        Data matrix containing points drawn from a union of subspaces as its rows
    label : shape (num_subspaces * num_points_per_subspace)
        Membership of each data point to the subspace it lies in
    """

    data = np.empty((num_points_per_subspace* num_subspaces, ambient_dim))
    label = np.empty(num_points_per_subspace * num_subspaces, dtype=int)
  
    for i in range(num_subspaces):
        basis = np.random.normal(size=(ambient_dim, subspace_dim))
        basis = orth(basis)
        coeff = np.random.normal(size=(subspace_dim, num_points_per_subspace))
        coeff = normalize(coeff, norm='l2', axis=0, copy=False)
        data_per_subspace = np.matmul(basis, coeff).T

        base_index = i*num_points_per_subspace
        data[(0+base_index):(num_points_per_subspace+base_index), :] = data_per_subspace
        label[0+base_index:num_points_per_subspace+base_index,] = i

    data += np.random.normal(size=(num_points_per_subspace * num_subspaces, ambient_dim)) * noise_level
  
    return data, label