biocore / gemelli / 127

pending completion

Build # 127

Build Type

Pull #36

travis-ci

Committed by

web-flow

Commit Message

<a href="https://github.com/biocore/gemelli/commit/<a class=hub.com/biocore/gemelli/commit/<a class="double-link" href="https://git"><a class=hub.com/biocore/gemelli/commit/<a class="double-link" href="https://git"><a class=hub.com/biocore/gemelli/commit/<a class="double-link" href="https://git"><a class=hub.com/biocore/gemelli/commit/62ec61b8db7e74de0b3abfbece8a6af38b83ee6e">62ec61b8d">&lt;a href=&quot;https://github.com/biocore/gemelli/commit/</a><a class="double-link" href="https://github.com/biocore/gemelli/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/biocore/gemelli/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/biocore/gemelli/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/biocore/gemelli/commit/&lt;a class=&quot;double-link&quot; href=&quot;https://git">&lt;a class=</a>hub.com/biocore/gemelli/commit/62ec61b8db7e74de0b3abfbece8a6af38b83ee6e">62ec61b8d</a><a href="https://github.com/biocore/gemelli/commit/62ec61b8db7e74de0b3abfbece8a6af38b83ee6e">&quot;&gt;&amp;lt;a href=&amp;quot;https://github.com/biocore/gemelli/commit/&lt;/a&gt;&lt;a class=&quot;double-link&quot; href=&quot;https://github.com/biocore/gemelli/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/biocore/gemelli/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/biocore/gemelli/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/biocore/gemelli/commit/&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://git&quot;&gt;&amp;lt;a class=&lt;/a&gt;hub.com/biocore/gemelli/commit/62ec61b8db7e74de0b3abfbece8a6af38b83ee6e&quot;&gt;62ec61b8d&lt;/a&gt;&lt;a href=&quot;https://github.com/biocore/gemelli/commit/62ec61b8db7e74de0b3abfbece8a6af38b83ee6e&quot;&gt;&amp;quot;&amp;gt;&amp;amp;lt;a href=&amp;amp;quot;https://github.com/biocore/gemelli/commit/&amp;lt;/a&amp;gt;&amp;lt;a class=&amp;quot;double-link&amp;quot; href=&amp;quot;https://github.com/biocore/gemelli/commit/&amp;amp;lt;a class=&amp;amp;quot;double-link&am... (continued)

Pull Request Pull Request #36: Add RPCA and standalone RCLR functionality

Run Details

338 of 384 branches covered (88.02%)

Branch coverage included in aggregate %.

344 of 344 new or added lines in 11 files covered. (100.0%)

860 of 905 relevant lines covered (95.03%)

0.95 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

92.86

/gemelli/optspace.py

# ----------------------------------------------------------------------------
# Copyright (c) 2019--, gemelli development team.
#
# Distributed under the terms of the Modified BSD License.
#
# The full license is in the file COPYING.txt, distributed with this software.
# ----------------------------------------------------------------------------

import warnings
import numpy as np
from numpy.matlib import repmat
from numpy.linalg import norm
from scipy.sparse.linalg import svds


class OptSpace(object):
    """

    OptSpace is a matrix completion algorithm based on a singular value
    decomposition (SVD) optimized on a local manifold. It has been shown to
    be quite robust to noise in low rank datasets (1).
    The objective function that it is trying to optimize
    over is given by:

        min(P|(Y-U*S*V^{T})|_{2}^{2}

    U and V are matrices that are trying to be estimated and S
    is analogous to a matrix of eigenvalues. Y are the
    observed values and P is a function such that
    the errors between Y and USV are only computed
    on the nonzero entries.

    """

    def __init__(
            self,
            n_components,
            max_iterations,
            tol,
            step_size=10000,
            resolution_limit=20,
            sign=-1):
        """

        Parameters
        ----------
        obs: numpy.ndarray - a rclr preprocessed matrix of shape (M,N)
            with missing values set to zero or np.nan.
            N = Features (i.e. OTUs, metabolites)
            M = Samples

        n_components: int or {"optspace"}, optional
            The underlying rank of the dataset.
            This will also control the number of components
            (axis) in the U and V loadings. The value can either
            be hard set by the user or estimated through the
            optspace algorithm.

        max_iterations: int
            The maximum number of convex iterations to optimize the solution
            If iteration is not specified, then the default iteration is 5.
            Which redcues to a satisfactory error threshold.

        tol: float
            Error reduction break, if the error reduced is
            less than this value it will return the solution

        step_size: int, optional : Default is 10000
            The gradient decent step size, this will be
            optimized by the line search.

        resolution_limit: int, optional : Default is 20
            The gradient decent line search resolution limit.
            Where the resolution is line > 2**-resolution_limit.

        sign: int, optional : Default is -1
            Can be one or negative one. This controls the
            sign correction in the gradient decent U, V
            updates.

        Returns
        -------
        self.U: numpy.ndarray - "Sample Loadings" or the unitary matrix
            having left singular vectors as columns.
            Of shape (M, n_components)

        self.s: numpy.ndarray - The singular values,
            sorted in non-increasing order.
            Of shape (n_components, n_components).

        self.V: numpy.ndarray - "Feature Loadings" or Unitary matrix
            having right singular vectors as rows.
            Of shape (N, n_components)

        References
        ----------
        .. [1] Keshavan RH, Oh S, Montanari A. 2009. Matrix completion
                from a few entries (2009_ IEEE International
                Symposium on Information Theory
        """

        self.n_components = n_components
        self.max_iterations = max_iterations
        self.tol = tol
        self.step_size = step_size
        self.resolution_limit = resolution_limit
        self.sign = sign

    def solve(self, obs):
        # adjust iteration indexing by one
        self.max_iterations += 1
        # Convert any nan input to zeros
        # optspace considers zero  and only zero
        # as missing.
        obs[np.isnan(obs)] = 0
        # generate a mask that tracks where missing
        # values exist in the obs dataset
        mask = (np.abs(obs) > 0).astype(np.int)
        # save the shape of the matrix
        n, m = obs.shape
        # get a measure of sparsity level
        total_nonzeros = np.count_nonzero(mask)
        eps = total_nonzeros / np.sqrt(m * n)
        if isinstance(self.n_components, str):
            if self.n_components.lower() == 'auto':
                # estimate the rank of the matrix
                self.n_components = rank_estimate(obs, eps)
                # check estimate again
                if self.n_components >= min(n, m) - 1:
                    warnings.warn('Your matrix is estimated '
                                  'to be high-rank.',
                                  RuntimeWarning)
                # print rank estimate
                print('Estimated rank is %i' % self.n_components)
            else:
                raise ValueError("n-components must be an "
                                 "integer or 'auto'.")
        # raise future warning if hard set
        elif isinstance(self.n_components, int):
            if self.n_components > (min(n, m) - 1):
                raise ValueError("n-components must be at most"
                                 " 1 minus the min. shape of the"
                                 " input matrix.")
        # otherwise rase an error.
        else:
            raise ValueError("n-components must be "
                             "an interger or 'auto'")
        # The rescaling factor compensates the smaller average size of
        # the of the missing entries (mask) with respect to obs
        rescal_param = np.count_nonzero(mask) * self.n_components
        rescal_param = np.sqrt(rescal_param / (norm(obs, 'fro') ** 2))
        obs = obs * rescal_param
        # Our initial first guess are the loadings generated
        # by the traditional SVD
        U, S, V = svds(obs, self.n_components, which='LM')
        # the shape and number of non-zero values
        # can set the input perameters for the gradient
        # decent
        rho = eps * n
        U = U * np.sqrt(n)
        V = (V * np.sqrt(m)).T
        S = S / eps
        # generate the new singular values from
        # the initialization of U and V
        S = singular_values(U, V, obs, mask)
        # initialize the difference between the
        # observed values of the matrix and the
        # the imputed matrix generated by the loadings
        # from this point on we call this "distortion"
        obs_error = obs - U.dot(S).dot(V.T)
        # initialize the distortion matrix of line above
        dist = np.zeros(self.max_iterations + 1)
        # starting initialization of the distortion between obs and imputed
        dist[0] = norm(np.multiply(obs_error, mask), 'fro') / \
            np.sqrt(total_nonzeros)
        # we will perform gradient decent for at most self.max_iterations
        for i in range(1, self.max_iterations):
            # generate new optimized loadings from F(X,Y) [1]
            U_update, V_update = gradient_decent(
                U, V, S, obs, mask, self.step_size, rho)
            # Line search for the optimum jump length in gradient decent
            line = line_search(
                U,
                U_update,
                V,
                V_update,
                S,
                obs,
                mask,
                self.step_size,
                rho,
                resolution_limit=self.resolution_limit)
            # with line and iterations loading update U,V loadings
            U = U - self.sign * line * U_update
            V = V - self.sign * line * V_update
            # generate new singular values from the new
            # loadings
            S = singular_values(U, V, obs, mask)
            # Compute the distortion
            obs_error = obs - U.dot(S).dot(V.T)
            # update the new distortion
            dist[i + 1] = norm(np.multiply(obs_error, mask),
                               'fro') / np.sqrt(total_nonzeros)
            # if the gradient decent has coverged then break the loop
            # and return the results
            if(dist[i + 1] < self.tol):
                break
        # compensates the smaller average size of
        # observed values vs. missing
        S = S / rescal_param
        # ensure the loadings are sorted properly
        U, S, V = svd_sort(U, S, V)
        return U, S, V


def svd_sort(U, S, V):
    """
    Sorting via the s matrix from SVD. In addition to
    sign correction from the U matrix to ensure a
    deterministic output.

    Parameters
    ----------
    U: array-like
        U matrix from SVD
    V: array-like
        V matrix from SVD
    S: array-like
        S matrix from SVD

    Notes
    -----
    S matrix can be off diagonal elements.
    """
    # See https://github.com/scikit-learn/scikit-learn/
    # blob/7b136e92acf49d46251479b75c88cba632de1937/sklearn/
    # decomposition/pca.py#L510-#L518 for context.
    # Because svds do not abide by the normal
    # conventions in scipy.linalg.svd/randomized_svd
    # the output has to be reversed
    idx = np.argsort(np.diag(S))[::-1]
    # sorting following the solution
    # provided by https://stackoverflow.com/
    # questions/36381356/sort-matrix-based
    # -on-its-diagonal-entries
    S = S[idx, :][:, idx]
    U = U[:, idx]
    V = V[:, idx]
    return U, S, V


def cost_function(U, V, S, obs, mask, step_size, rho):
    """
    Parameters
    ----------
    U, V, S, obs, mask, step_size, rho

    Notes
    -----
    M ~ imputed
    """
    # shape of subject loading
    n, n_components = U.shape
    # calculate the Frobenius norm between observed values and imputed values
    distortion = np.sum(
        np.sum((np.multiply((U.dot(S).dot(V.T) - obs), mask)**2))) / 2
    # calculate the  cartesian product of two Grassmann manifolds
    V_manifold = rho * grassmann_manifold_one(V, step_size, n_components)
    U_manifold = rho * grassmann_manifold_one(U, step_size, n_components)
    return distortion + V_manifold + U_manifold


def gradient_decent(U, V, S, obs, mask, step_size, rho):
    """
    A single iteration of the gradient decent update to the U and V matrix.

    Parameters
    ----------
    U, V, S, obs, mask, step_size, rho

    Returns
    -------
    U_update, V_update
    """
    # shape of loadings
    n, n_components = U.shape
    m, n_components = V.shape
    # generate the new values imputed
    # from the loadings from obs values
    US = U.dot(S)
    VS = V.dot(S.T)
    imputed = US.dot(V.T)
    # calculate the U and V distortion
    Qu = U.T.dot(np.multiply((obs - imputed), mask)).dot(VS) / n
    Qv = V.T.dot(np.multiply((obs - imputed), mask).T).dot(US) / m
    # create new loadings based on the distortion between obs and imputed
    # pass these loadings to back to the decent.
    U_update = np.multiply((imputed - obs), mask).dot(VS) + U.dot(Qu) + \
        rho * grassmann_manifold_two(U, step_size, n_components)
    V_update = np.multiply((imputed - obs), mask).T.dot(US) + V.dot(Qv) + \
        rho * grassmann_manifold_two(V, step_size, n_components)
    return U_update, V_update


def line_search(
        U,
        U_update,
        V,
        V_update,
        S,
        obs,
        mask,
        step_size,
        rho,
        resolution_limit=20,
        line=-1e-1):
    """
    An exact line search
        gradient decent converging for a quadratic function

    Parameters
    ----------
    U, U_update, V, V_update, S, obs, mask, step_size, rho
    """

    norm_update = norm(U_update, 'fro')**2 + norm(V_update, 'fro')**2
    # this is the resolution limit (line > 2**-20)
    cost = np.zeros(resolution_limit + 1)
    cost[0] = cost_function(U, V, S, obs, mask, step_size, rho)
    for i in range(resolution_limit):
        cost[i +
             1] = cost_function(U +
                                line *
                                U_update, V +
                                line *
                                V_update, S, obs, mask, step_size, rho)
        if((cost[i + 1] - cost[0]) <= .5 * line * norm_update):
            return line
        line = line / 2
    return line


def singular_values(U, V, obs, mask):
    """
    Generates the singular values from the updated
    U & V loadings for one iteration.

    Parameters
    ----------
    U, V, obs, mask

    """

    n, n_components = U.shape
    C = np.ravel(U.T.dot(obs).dot(V))
    A = np.zeros((n_components * n_components, n_components * n_components))
    for i in range(n_components):
        for j in range(n_components):
            ind = j * n_components + i
            x = U[:, i].reshape(1, len(U[:, i]))
            manifold = V[:, j].reshape(len(V[:, j]), 1)
            tmp = np.multiply((manifold.dot(x)).T, mask)
            temp = U.T.dot(tmp).dot(V)
            A[:, ind] = np.ravel(temp)
    S = np.linalg.lstsq(A, C, rcond=1e-12)[0]
    S = S.reshape((n_components, n_components)).T

    return S


def grassmann_manifold_one(U, step_size, n_components):
    """
    The first Grassmann Manifold

    Parameters
    ----------
    U, step_size, n_components
    """
    # get the step from the manifold
    step = np.sum(U**2, axis=1) / (2 * step_size * n_components)
    manifold = np.exp((step - 1)**2) - 1
    manifold[step < 1] = 0
    manifold[manifold == np.inf] = 0
    return manifold.sum()


def grassmann_manifold_two(U, step_size, n_components):
    """
    The second Grassmann Manifold

    Parameters
    ----------
    U, step_size, n_components


    """
    # get the step from the manifold
    step = np.sum(U**2, axis=1) / (2 * step_size * n_components)
    step = 2 * np.multiply(np.exp((step - 1)**2), (step - 1))
    step[step < 0] = 0
    step = step.reshape(len(step), 1)
    step = np.multiply(U, repmat(step, 1, n_components)) / \
        (step_size * n_components)
    return step


def rank_estimate(obs, eps, k=20, lam=0.05,
                  min_rank=3, max_iter=5000):

    """
    This function estimates the rank of a
    sparse matrix (i.e. with missing values).

    Parameters
    ----------
    obs: numpy.ndarray - a rclr preprocessed matrix of shape (M,N)
        with missing values set to zero or np.nan.
        N = Features (i.e. OTUs, metabolites)
        M = Samples

    eps: float - Measure of the level of sparsity
        Equivalent to obs N-non-zeros / sqrt(obs.shape)

    k: int - Max number of singular values / rank

    lam: float - Step in the iteration

    min_rank: int - The min. rank allowed

    Returns
    -------
    The estimated rank of the matrix.

    References
    ----------
    .. [1] Part C in Keshavan, R. H., Montanari,
           A. & Oh, S. Low-rank matrix completion
           with noisy observations: A quantitative
           comparison. in 2009 47th Annual Allerton
           Conference on Communication, Control,
           and Computing (Allerton) 1216–1222 (2009).
    """

    # dim. of the data
    n, m = obs.shape
    # get N-singular values
    s = svds(obs,  min(k, n, m) - 1, which='LM',
             return_singular_vectors=False)[::-1]
    # get N+1 singular values
    s_one = s[:-1] - s[1:]
    # simplify iterations
    s_one_ = s_one / np.mean(s_one[-10:])
    # iteration one
    r_one = 0
    iter_ = 0
    while r_one <= 0:
        cost = []
        # get the cost
        for idx in range(s_one_.shape[0]):
            cost.append(lam * max(s_one_[idx:]) + idx)
        # estimate the rank
        r_one = np.argmin(cost)
        lam += lam
        iter_ += 1
        if iter_ > max_iter:
            break

    # iteration two
    cost = []
    # estimate the rank
    for idx in range(s.shape[0] - 1):
        cost.append((s[idx + 1] + np.sqrt(idx * eps) * s[0] / eps) / s[idx])
    r_two = np.argmin(cost)
    # return the final estimate
    return max(r_one, r_two, min_rank)

1	# ----------------------------------------------------------------------------
2	# Copyright (c) 2019--, gemelli development team.
3	#
4	# Distributed under the terms of the Modified BSD License.
5	#
6	# The full license is in the file COPYING.txt, distributed with this software.
7	# ----------------------------------------------------------------------------
8
9	import warnings	1✔
10	import numpy as np	1✔
11	from numpy.matlib import repmat	1✔
12	from numpy.linalg import norm	1✔
13	from scipy.sparse.linalg import svds	1✔
14
15
16	class OptSpace(object):	1✔
17	"""
18
19	OptSpace is a matrix completion algorithm based on a singular value
20	decomposition (SVD) optimized on a local manifold. It has been shown to
21	be quite robust to noise in low rank datasets (1).
22	The objective function that it is trying to optimize
23	over is given by:
24
25	min(P\|(Y-USV^{T})\|_{2}^{2}
26
27	U and V are matrices that are trying to be estimated and S
28	is analogous to a matrix of eigenvalues. Y are the
29	observed values and P is a function such that
30	the errors between Y and USV are only computed
31	on the nonzero entries.
32
33	"""
34
35	def __init__(	1✔
36	self,
37	n_components,
38	max_iterations,
39	tol,
40	step_size=10000,
41	resolution_limit=20,
42	sign=-1):
43	"""
44
45	Parameters
46	----------
47	obs: numpy.ndarray - a rclr preprocessed matrix of shape (M,N)
48	with missing values set to zero or np.nan.
49	N = Features (i.e. OTUs, metabolites)
50	M = Samples
51
52	n_components: int or {"optspace"}, optional
53	The underlying rank of the dataset.
54	This will also control the number of components
55	(axis) in the U and V loadings. The value can either
56	be hard set by the user or estimated through the
57	optspace algorithm.
58
59	max_iterations: int
60	The maximum number of convex iterations to optimize the solution
61	If iteration is not specified, then the default iteration is 5.
62	Which redcues to a satisfactory error threshold.
63
64	tol: float
65	Error reduction break, if the error reduced is
66	less than this value it will return the solution
67
68	step_size: int, optional : Default is 10000
69	The gradient decent step size, this will be
70	optimized by the line search.
71
72	resolution_limit: int, optional : Default is 20
73	The gradient decent line search resolution limit.
74	Where the resolution is line > 2**-resolution_limit.
75
76	sign: int, optional : Default is -1
77	Can be one or negative one. This controls the
78	sign correction in the gradient decent U, V
79	updates.
80
81	Returns
82	-------
83	self.U: numpy.ndarray - "Sample Loadings" or the unitary matrix
84	having left singular vectors as columns.
85	Of shape (M, n_components)
86
87	self.s: numpy.ndarray - The singular values,
88	sorted in non-increasing order.
89	Of shape (n_components, n_components).
90
91	self.V: numpy.ndarray - "Feature Loadings" or Unitary matrix
92	having right singular vectors as rows.
93	Of shape (N, n_components)
94
95	References
96	----------
97	.. [1] Keshavan RH, Oh S, Montanari A. 2009. Matrix completion
98	from a few entries (2009_ IEEE International
99	Symposium on Information Theory
100	"""
101
102	self.n_components = n_components	1✔
103	self.max_iterations = max_iterations	1✔
104	self.tol = tol	1✔
105	self.step_size = step_size	1✔
106	self.resolution_limit = resolution_limit	1✔
107	self.sign = sign	1✔
108
109	def solve(self, obs):	1✔
110	# adjust iteration indexing by one
111	self.max_iterations += 1	1✔
112	# Convert any nan input to zeros
113	# optspace considers zero and only zero
114	# as missing.
115	obs[np.isnan(obs)] = 0	1✔
116	# generate a mask that tracks where missing
117	# values exist in the obs dataset
118	mask = (np.abs(obs) > 0).astype(np.int)	1✔
119	# save the shape of the matrix
120	n, m = obs.shape	1✔
121	# get a measure of sparsity level
122	total_nonzeros = np.count_nonzero(mask)	1✔
123	eps = total_nonzeros / np.sqrt(m * n)	1✔
124	if isinstance(self.n_components, str):	1✔
125	if self.n_components.lower() == 'auto':	1!
126	# estimate the rank of the matrix
127	self.n_components = rank_estimate(obs, eps)	1✔
128	# check estimate again
129	if self.n_components >= min(n, m) - 1:	1!
130	warnings.warn('Your matrix is estimated '	×
131	'to be high-rank.',
132	RuntimeWarning)
133	# print rank estimate
134	print('Estimated rank is %i' % self.n_components)	1✔
135	else:
136	raise ValueError("n-components must be an "	×
137	"integer or 'auto'.")
138	# raise future warning if hard set
139	elif isinstance(self.n_components, int):	1!
140	if self.n_components > (min(n, m) - 1):	1!
141	raise ValueError("n-components must be at most"	×
142	" 1 minus the min. shape of the"
143	" input matrix.")
144	# otherwise rase an error.
145	else:
146	raise ValueError("n-components must be "	×
147	"an interger or 'auto'")
148	# The rescaling factor compensates the smaller average size of
149	# the of the missing entries (mask) with respect to obs
150	rescal_param = np.count_nonzero(mask) * self.n_components	1✔
151	rescal_param = np.sqrt(rescal_param / (norm(obs, 'fro') ** 2))	1✔
152	obs = obs * rescal_param	1✔
153	# Our initial first guess are the loadings generated
154	# by the traditional SVD
155	U, S, V = svds(obs, self.n_components, which='LM')	1✔
156	# the shape and number of non-zero values
157	# can set the input perameters for the gradient
158	# decent
159	rho = eps * n	1✔
160	U = U * np.sqrt(n)	1✔
161	V = (V * np.sqrt(m)).T	1✔
162	S = S / eps	1✔
163	# generate the new singular values from
164	# the initialization of U and V
165	S = singular_values(U, V, obs, mask)	1✔
166	# initialize the difference between the
167	# observed values of the matrix and the
168	# the imputed matrix generated by the loadings
169	# from this point on we call this "distortion"
170	obs_error = obs - U.dot(S).dot(V.T)	1✔
171	# initialize the distortion matrix of line above
172	dist = np.zeros(self.max_iterations + 1)	1✔
173	# starting initialization of the distortion between obs and imputed
174	dist[0] = norm(np.multiply(obs_error, mask), 'fro') / \	1✔
175	np.sqrt(total_nonzeros)
176	# we will perform gradient decent for at most self.max_iterations
177	for i in range(1, self.max_iterations):	1✔
178	# generate new optimized loadings from F(X,Y) [1]
179	U_update, V_update = gradient_decent(	1✔
180	U, V, S, obs, mask, self.step_size, rho)
181	# Line search for the optimum jump length in gradient decent
182	line = line_search(	1✔
183	U,
184	U_update,
185	V,
186	V_update,
187	S,
188	obs,
189	mask,
190	self.step_size,
191	rho,
192	resolution_limit=self.resolution_limit)
193	# with line and iterations loading update U,V loadings
194	U = U - self.sign * line * U_update	1✔
195	V = V - self.sign * line * V_update	1✔
196	# generate new singular values from the new
197	# loadings
198	S = singular_values(U, V, obs, mask)	1✔
199	# Compute the distortion
200	obs_error = obs - U.dot(S).dot(V.T)	1✔
201	# update the new distortion
202	dist[i + 1] = norm(np.multiply(obs_error, mask),	1✔
203	'fro') / np.sqrt(total_nonzeros)
204	# if the gradient decent has coverged then break the loop
205	# and return the results
206	if(dist[i + 1] < self.tol):	1!
207	break	×
208	# compensates the smaller average size of
209	# observed values vs. missing
210	S = S / rescal_param	1✔
211	# ensure the loadings are sorted properly
212	U, S, V = svd_sort(U, S, V)	1✔
213	return U, S, V	1✔
214
215
216	def svd_sort(U, S, V):	1✔
217	"""
218	Sorting via the s matrix from SVD. In addition to
219	sign correction from the U matrix to ensure a
220	deterministic output.
221
222	Parameters
223	----------
224	U: array-like
225	U matrix from SVD
226	V: array-like
227	V matrix from SVD
228	S: array-like
229	S matrix from SVD
230
231	Notes
232	-----
233	S matrix can be off diagonal elements.
234	"""
235	# See https://github.com/scikit-learn/scikit-learn/
236	# blob/7b136e92acf49d46251479b75c88cba632de1937/sklearn/
237	# decomposition/pca.py#L510-#L518 for context.
238	# Because svds do not abide by the normal
239	# conventions in scipy.linalg.svd/randomized_svd
240	# the output has to be reversed
241	idx = np.argsort(np.diag(S))[::-1]	1✔
242	# sorting following the solution
243	# provided by https://stackoverflow.com/
244	# questions/36381356/sort-matrix-based
245	# -on-its-diagonal-entries
246	S = S[idx, :][:, idx]	1✔
247	U = U[:, idx]	1✔
248	V = V[:, idx]	1✔
249	return U, S, V	1✔
250
251
252	def cost_function(U, V, S, obs, mask, step_size, rho):	1✔
253	"""
254	Parameters
255	----------
256	U, V, S, obs, mask, step_size, rho
257
258	Notes
259	-----
260	M ~ imputed
261	"""
262	# shape of subject loading
263	n, n_components = U.shape	1✔
264	# calculate the Frobenius norm between observed values and imputed values
265	distortion = np.sum(	1✔
266	np.sum((np.multiply((U.dot(S).dot(V.T) - obs), mask)**2))) / 2
267	# calculate the cartesian product of two Grassmann manifolds
268	V_manifold = rho * grassmann_manifold_one(V, step_size, n_components)	1✔
269	U_manifold = rho * grassmann_manifold_one(U, step_size, n_components)	1✔
270	return distortion + V_manifold + U_manifold	1✔
271
272
273	def gradient_decent(U, V, S, obs, mask, step_size, rho):	1✔
274	"""
275	A single iteration of the gradient decent update to the U and V matrix.
276
277	Parameters
278	----------
279	U, V, S, obs, mask, step_size, rho
280
281	Returns
282	-------
283	U_update, V_update
284	"""
285	# shape of loadings
286	n, n_components = U.shape	1✔
287	m, n_components = V.shape	1✔
288	# generate the new values imputed
289	# from the loadings from obs values
290	US = U.dot(S)	1✔
291	VS = V.dot(S.T)	1✔
292	imputed = US.dot(V.T)	1✔
293	# calculate the U and V distortion
294	Qu = U.T.dot(np.multiply((obs - imputed), mask)).dot(VS) / n	1✔
295	Qv = V.T.dot(np.multiply((obs - imputed), mask).T).dot(US) / m	1✔
296	# create new loadings based on the distortion between obs and imputed
297	# pass these loadings to back to the decent.
298	U_update = np.multiply((imputed - obs), mask).dot(VS) + U.dot(Qu) + \	1✔
299	rho * grassmann_manifold_two(U, step_size, n_components)
300	V_update = np.multiply((imputed - obs), mask).T.dot(US) + V.dot(Qv) + \	1✔
301	rho * grassmann_manifold_two(V, step_size, n_components)
302	return U_update, V_update	1✔
303
304
305	def line_search(	1✔
306	U,
307	U_update,
308	V,
309	V_update,
310	S,
311	obs,
312	mask,
313	step_size,
314	rho,
315	resolution_limit=20,
316	line=-1e-1):
317	"""
318	An exact line search
319	gradient decent converging for a quadratic function
320
321	Parameters
322	----------
323	U, U_update, V, V_update, S, obs, mask, step_size, rho
324	"""
325
326	norm_update = norm(U_update, 'fro')2 + norm(V_update, 'fro')2	1✔
327	# this is the resolution limit (line > 2**-20)
328	cost = np.zeros(resolution_limit + 1)	1✔
329	cost[0] = cost_function(U, V, S, obs, mask, step_size, rho)	1✔
330	for i in range(resolution_limit):	1✔
331	cost[i +	1✔
332	1] = cost_function(U +
333	line *
334	U_update, V +
335	line *
336	V_update, S, obs, mask, step_size, rho)
337	if((cost[i + 1] - cost[0]) <= .5 * line * norm_update):	1✔
338	return line	1✔
339	line = line / 2	1✔
340	return line	1✔
341
342
343	def singular_values(U, V, obs, mask):	1✔
344	"""
345	Generates the singular values from the updated
346	U & V loadings for one iteration.
347
348	Parameters
349	----------
350	U, V, obs, mask
351
352	"""
353
354	n, n_components = U.shape	1✔
355	C = np.ravel(U.T.dot(obs).dot(V))	1✔
356	A = np.zeros((n_components * n_components, n_components * n_components))	1✔
357	for i in range(n_components):	1✔
358	for j in range(n_components):	1✔
359	ind = j * n_components + i	1✔
360	x = U[:, i].reshape(1, len(U[:, i]))	1✔
361	manifold = V[:, j].reshape(len(V[:, j]), 1)	1✔
362	tmp = np.multiply((manifold.dot(x)).T, mask)	1✔
363	temp = U.T.dot(tmp).dot(V)	1✔
364	A[:, ind] = np.ravel(temp)	1✔
365	S = np.linalg.lstsq(A, C, rcond=1e-12)[0]	1✔
366	S = S.reshape((n_components, n_components)).T	1✔
367
368	return S	1✔
369
370
371	def grassmann_manifold_one(U, step_size, n_components):	1✔
372	"""
373	The first Grassmann Manifold
374
375	Parameters
376	----------
377	U, step_size, n_components
378	"""
379	# get the step from the manifold
380	step = np.sum(U*2, axis=1) / (2 step_size * n_components)	1✔
381	manifold = np.exp((step - 1)**2) - 1	1✔
382	manifold[step < 1] = 0	1✔
383	manifold[manifold == np.inf] = 0	1✔
384	return manifold.sum()	1✔
385
386
387	def grassmann_manifold_two(U, step_size, n_components):	1✔
388	"""
389	The second Grassmann Manifold
390
391	Parameters
392	----------
393	U, step_size, n_components
394
395
396	"""
397	# get the step from the manifold
398	step = np.sum(U*2, axis=1) / (2 step_size * n_components)	1✔
399	step = 2 * np.multiply(np.exp((step - 1)**2), (step - 1))	1✔
400	step[step < 0] = 0	1✔
401	step = step.reshape(len(step), 1)	1✔
402	step = np.multiply(U, repmat(step, 1, n_components)) / \	1✔
403	(step_size * n_components)
404	return step	1✔
405
406
407	def rank_estimate(obs, eps, k=20, lam=0.05,	1✔
408	min_rank=3, max_iter=5000):
409
410	"""
411	This function estimates the rank of a
412	sparse matrix (i.e. with missing values).
413
414	Parameters
415	----------
416	obs: numpy.ndarray - a rclr preprocessed matrix of shape (M,N)
417	with missing values set to zero or np.nan.
418	N = Features (i.e. OTUs, metabolites)
419	M = Samples
420
421	eps: float - Measure of the level of sparsity
422	Equivalent to obs N-non-zeros / sqrt(obs.shape)
423
424	k: int - Max number of singular values / rank
425
426	lam: float - Step in the iteration
427
428	min_rank: int - The min. rank allowed
429
430	Returns
431	-------
432	The estimated rank of the matrix.
433
434	References
435	----------
436	.. [1] Part C in Keshavan, R. H., Montanari,
437	A. & Oh, S. Low-rank matrix completion
438	with noisy observations: A quantitative
439	comparison. in 2009 47th Annual Allerton
440	Conference on Communication, Control,
441	and Computing (Allerton) 1216–1222 (2009).
442	"""
443
444	# dim. of the data
445	n, m = obs.shape	1✔
446	# get N-singular values
447	s = svds(obs, min(k, n, m) - 1, which='LM',	1✔
448	return_singular_vectors=False)[::-1]
449	# get N+1 singular values
450	s_one = s[:-1] - s[1:]	1✔
451	# simplify iterations
452	s_one_ = s_one / np.mean(s_one[-10:])	1✔
453	# iteration one
454	r_one = 0	1✔
455	iter_ = 0	1✔
456	while r_one <= 0:	1✔
457	cost = []	1✔
458	# get the cost
459	for idx in range(s_one_.shape[0]):	1✔
460	cost.append(lam * max(s_one_[idx:]) + idx)	1✔
461	# estimate the rank
462	r_one = np.argmin(cost)	1✔
463	lam += lam	1✔
464	iter_ += 1	1✔
465	if iter_ > max_iter:	1!
466	break	×
467
468	# iteration two
469	cost = []	1✔
470	# estimate the rank
471	for idx in range(s.shape[0] - 1):	1✔
472	cost.append((s[idx + 1] + np.sqrt(idx * eps) * s[0] / eps) / s[idx])	1✔
473	r_two = np.argmin(cost)	1✔
474	# return the final estimate
475	return max(r_one, r_two, min_rank)	1✔

biocore / gemelli / 127

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous