8000 Merge pull request #139 from lan496/reproducibility-for-tucker · tensorly/tensorly@dc67ac6 · GitHub
[go: up one dir, main page]
Skip to content

Commit dc67ac6

Browse files
JeanKossaifiJeanKossaifi
authored
Merge pull request #139 from lan496/reproducibility-for-tucker
Making partial_svd and tucker reproducible
2 parents f1a7bf2 + 2a38abf commit dc67ac6

File tree

9 files changed

+79
-25
lines changed
  • tests
  • tests

    • 9 files changed

    +79
    -25
    lines changed

    ‎tensorly/backend/core.py

    Lines changed: 24 additions & 8 deletions
    Original file line numberDiff line numberDiff line change
    @@ -552,10 +552,10 @@ def stack(arrays, axis=0):
    552552
    Join a sequence of arrays along a new axis.
    553553
    """
    554554
    raise NotImplementedError
    555-
    555+
    556556
    def eps(self, dtype):
    557557
    return self.finfo(dtype).eps
    558-
    558+
    559559
    def finfo(self, dtype):
    560560
    return np.finfo(self.to_numpy(self.tensor([], dtype=dtype)).dtype)
    561561

    @@ -632,12 +632,12 @@ def kr(self, matrices, weights=None, mask=None):
    632632
    a = self.reshape(res, (s1, 1, s2))
    633633
    b = self.reshape(e, (1, s3, s4))
    634634
    res = self.reshape(a * b, (-1, n_col))
    635-
    635+
    636636
    m = self.reshape(mask, (-1, 1)) if mask is not None else 1
    637-
    637+
    638638
    return res*m
    639639

    640-
    def partial_svd(self, matrix, n_eigenvecs=None):
    640+
    def partial_svd(self, matrix, n_eigenvecs=None, random_state=None, **kwargs):
    641641
    """Computes a fast partial SVD on `matrix`
    642642
    643643
    If `n_eigenvecs` is specified, sparse eigendecomposition is used on
    @@ -649,6 +649,10 @@ def partial_svd(self, matrix, n_eigenvecs=None):
    649649
    A 2D tensor.
    650650
    n_eigenvecs : int, optional, default is None
    651651
    If specified, number of eigen[vectors-values] to return.
    652+
    random_state: {None, int, np.random.RandomState}
    653+
    If specified, use it for sampling starting vector in a partial SVD(scipy.sparse.linalg.eigsh)
    654+
    **kwargs : optional
    655+
    kwargs are used to absorb the difference of parameters among the other SVD functions
    652656
    653657
    Returns
    654658
    -------
    @@ -696,16 +700,28 @@ def partial_svd(self, matrix, n_eigenvecs=None):
    696700
    U, S, V = U[:, :n_eigenvecs], S[:n_eigenvecs], V[:n_eigenvecs, :]
    697701
    else:
    698702
    # We can perform a partial SVD
    703+
    # construct np.random.RandomState for sampling a starting vector
    704+
    if random_state is None:
    705+
    # if random_state is not specified, do not initialize a starting vector
    706+
    v0 = None
    707+
    elif isinstance(random_state, int):
    708+
    rns = np.random.RandomState(random_state)
    709+
    # initilize with [-1, 1] as in ARPACK
    710+
    v0 = rns.uniform(-1, 1, min_dim)
    711+
    elif isinstance(random_state, np.random.RandomState):
    712+
    # initilize with [-1, 1] as in ARPACK
    713+
    v0 = random_state.uniform(-1, 1, min_dim)
    714+
    699715
    # First choose whether to use X * X.T or X.T *X
    700716
    if dim_1 < dim_2:
    701717
    S, U = scipy.sparse.linalg.eigsh(
    702-
    np.dot(matrix, matrix.T.conj()), k=n_eigenvecs, which='LM'
    718+
    np.dot(matrix, matrix.T.conj()), k=n_eigenvecs, which='LM', v0=v0
    703719
    )
    704720
    S = np.where(np.abs(S) <= np.finfo(S.dtype).eps, 0, np.sqrt(S))
    705721
    V = np.dot(matrix.T.conj(), U * np.where(np.abs(S) <= np.finfo(S.dtype).eps, 0, 1/S)[None, :])
    706722
    else:
    707723
    S, V = scipy.sparse.linalg.eigsh(
    708-
    np.dot(matrix.T.conj(), matrix), k=n_eigenvecs, which='LM'
    724+
    np.dot(matrix.T.conj(), matrix), k=n_eigenvecs, which='LM', v0=v0
    709725
    )
    710726
    S = np.where(np.abs(S) <= np.finfo(S.dtype).eps, 0, np.sqrt(S))
    711727
    U = np.dot(matrix, V) * np.where(np.abs(S) <= np.finfo(S.dtype).eps, 0, 1/S)[None, :]
    @@ -718,4 +734,4 @@ def partial_svd(self, matrix, n_eigenvecs=None):
    718734
    U = self.tensor(U, **ctx)
    719735
    S = self.tensor(S, **ctx)
    720736
    V = self.tensor(V, **ctx)
    721-
    return U, S, V
    737+
    return U, S, V

    ‎tensorly/backend/cupy_backend.py

    Lines changed: 3 additions & 1 deletion
    Original file line numberDiff line numberDiff line change
    @@ -74,14 +74,16 @@ def solve(self, matrix1, matrix2):
    7474
    return self.tensor(res, **ctx)
    7575

    7676
    @staticmethod
    77-
    def truncated_svd(matrix, n_eigenvecs=None):
    77+
    def truncated_svd(matrix, n_eigenvecs=None, **kwargs):
    7878
    """Computes a truncated SVD on `matrix`
    7979
    8080
    Parameters
    8181
    ----------
    8282
    matrix : 2D-array
    8383
    n_eigenvecs : int, optional, default is None
    8484
    if specified, number of eigen[vectors-values] to return
    85+
    **kwargs : optional
    86+
    kwargs are used to absorb the difference of parameters among the other SVD functions
    8587
    8688
    Returns
    8789
    -------

    ‎tensorly/backend/mxnet_backend.py

    Lines changed: 3 additions & 1 deletion
    Original file line numberDiff line numberDiff line change
    @@ -223,7 +223,7 @@ def concatenate(tensors, axis):
    223223
    def stack(arrays, axis=0):
    224224
    return stack(*arrays, axis=axis)
    225225

    226-
    def symeig_svd(self, matrix, n_eigenvecs=None):
    226+
    def symeig_svd(self, matrix, n_eigenvecs=None, **kwargs):
    227227
    """Computes a truncated SVD on `matrix` using symeig
    228228
    229229
    Uses symeig on matrix.T.dot(matrix) or its transpose
    @@ -233,6 +233,8 @@ def symeig_svd(self, matrix, n_eigenvecs=None):
    233233
    matrix : 2D-array
    234234
    n_eigenvecs : int, optional, default is None
    235235
    if specified, number of eigen[vectors-values] to return
    236+
    **kwargs : optional
    237+
    kwargs are used to absorb the difference of parameters among the other SVD functions
    236238
    237239
    Returns
    238240
    -------

    ‎tensorly/backend/pytorch_backend.py

    Lines changed: 8 additions & 4 deletions
    Original file line numberDiff line numberDiff line change
    @@ -109,7 +109,7 @@ def solve(self, matrix1, matrix2):
    109109
    Notes
    110110
    -----
    111111
    112-
    Previously, this was implemented as follows::
    112+
    Previously, this was implemented as follows::
    113113
    114114
    if self.ndim(matrix2) < 2:
    115115
    # Currently, gesv doesn't support vectors for matrix2
    @@ -173,7 +173,7 @@ def argmax(input, axis=None):
    173173
    @staticmethod
    174174
    def stack(arrays, axis=0):
    175175
    return torch.stack(arrays, dim=axis)
    176-
    176+
    177177
    @staticmethod
    178178
    def conj(x, *args, **kwargs):
    179179
    """WARNING: IDENTITY FUNCTION (does nothing)
    @@ -201,14 +201,16 @@ def _reverse(tensor, axis=0):
    201201
    return tensor.index_select(axis, indices)
    202202

    203203
    @staticmethod
    204-
    def truncated_svd(matrix, n_eigenvecs=None):
    204+
    def truncated_svd(matrix, n_eigenvecs=None, **kwargs):
    205205
    """Computes a truncated SVD on `matrix` using pytorch's SVD
    206206
    207207
    Parameters
    208208
    ----------
    209209
    matrix : 2D-array
    210210
    n_eigenvecs : int, optional, default is None
    211211
    if specified, number of eigen[vectors-values] to return
    212+
    **kwargs : optional
    213+
    kwargs are used to absorb the difference of parameters among the other SVD functions
    212214
    213215
    111C Returns
    214216
    -------
    @@ -237,7 +239,7 @@ def truncated_svd(matrix, n_eigenvecs=None):
    237< F438 /td>239
    U, S, V = U[:, :n_eigenvecs], S[:n_eigenvecs], V[:n_eigenvecs, :]
    238240
    return U, S, V
    239241

    240-
    def symeig_svd(self, matrix, n_eigenvecs=None):
    242+
    def symeig_svd(self, matrix, n_eigenvecs=None, **kwargs):
    241243
    """Computes a truncated SVD on `matrix` using symeig
    242244
    243245
    Uses symeig on matrix.T.dot(matrix) or its transpose
    @@ -247,6 +249,8 @@ def symeig_svd(self, matrix, n_eigenvecs=None):
    247249
    matrix : 2D-array
    248250
    n_eigenvecs : int, optional, default is None
    249251
    if specified, number of eigen[vectors-values] to return
    252+
    **kwargs : optional
    253+
    kwargs are used to absorb the difference of parameters among the other SVD functions
    250254
    251255
    Returns
    252256
    -------

    ‎tensorly/backend/tensorflow_backend.py

    Lines changed: 4 additions & 2 deletions
    Original file line numberDiff line numberDiff line change
    @@ -98,7 +98,7 @@ def norm(tensor, order=2, axis=None):
    9898

    9999
    def dot(self, tensor1, tensor2):
    100100
    return tf.tensordot(tensor1, tensor2, axes=([self.ndim(tensor1) - 1], [0]))
    101-
    101+
    102102
    @staticmethod
    103103
    def conj(x, *args, **kwargs):
    104104
    """WARNING: IDENTITY FUNCTION (does nothing)
    @@ -119,14 +119,16 @@ def solve(lhs, rhs):
    119119
    return res
    120120

    121121
    @staticmethod
    122-
    def truncated_svd(matrix, n_eigenvecs=None):
    122+
    def truncated_svd(matrix, n_eigenvecs=None, **kwargs):
    123123
    """Computes an SVD on `matrix`
    124124
    125125
    Parameters
    126126
    ----------
    127127
    matrix : 2D-array
    128128
    n_eigenvecs : int, optional, default is None
    129129
    if specified, number of eigen[vectors-values] to return
    130+
    **kwargs : optional
    131+
    kwargs are used to absorb the difference of parameters among the other SVD functions
    130132
    131133
    Returns
    132134
    -------

    ‎tensorly/contrib/sparse/backend/numpy_backend.py

    Lines changed: 17 additions & 5 deletions
    Original file line numberDiff line numberDiff line change
    @@ -25,8 +25,8 @@ class NumpySparseBackend(Backend):
    2525
    backend_name = 'numpy.sparse'
    2626

    2727
    # moveaxis and shape are temporarily redefine to fix issue #131
    28-
    # Using the builting functionsn raises a TypeError:
    29-
    # no implementation found for 'numpy.shape' on types
    28+
    # Using the builting functionsn raises a TypeError:
    29+
    # no implementation found for 'numpy.shape' on types
    3030
    # that implement __array_function__: [<class 'sparse._coo.core.COO'>]
    3131
    def moveaxis(self, tensor, source, target):
    3232
    axes = list(range(self.ndim(tensor)))
    @@ -113,7 +113,7 @@ def solve(self, A, b):
    113113
    return x
    114114

    115115
    @staticmethod
    116-
    def partial_svd(matrix, n_eigenvecs=None):
    116+
    def partial_svd(matrix, n_eigenvecs=None, random_state=None, **kwargs):
    117117
    # Check that matrix is... a matrix!
    118118
    if matrix.ndim != 2:
    119119
    raise ValueError('matrix be a matrix. matrix.ndim is {} != 2'.format(
    @@ -150,6 +150,18 @@ def partial_svd(matrix, n_eigenvecs=None):
    150150
    if np.issubdtype(matrix.dtype, np.complexfloating):
    151151
    raise NotImplementedError("Complex dtypes")
    152152
    # We can perform a partial SVD
    153+
    # construct np.random.RandomState for sampling a starting vector
    154+
    if random_state is None:
    155+
    # if random_state is not specified, do not initialize a starting vector
    156+
    v0 = None
    157+
    elif isinstance(random_state, int):
    158+
    rns = np.random.RandomState(random_state)
    159+
    # initilize with [-1, 1] as in ARPACK
    160+
    v0 = rns.uniform(-1, 1, min_dim)
    161+
    elif isinstance(random_state, np.random.RandomState):
    162+
    # initilize with [-1, 1] as in ARPACK
    10000 163+
    v0 = random_state.uniform(-1, 1, min_dim)
    164+
    153165
    # First choose whether to use X * X.T or X.T *X
    154166
    if dim_1 < dim_2:
    155167
    conj = matrix.T
    @@ -160,7 +172,7 @@ def partial_svd(matrix, n_eigenvecs=None):
    160172
    # use dense form when sparse form will fail
    161173
    S, U = scipy.linalg.eigh(xxT.toarray())
    162174
    else:
    163-
    S, U = scipy.sparse.linalg.eigsh(xxT, k=n_eigenvecs, which='LM')
    175+
    S, U = scipy.sparse.linalg.eigsh(xxT, k=n_eigenvecs, which='LM', v0=v0)
    164176
    S = np.sqrt(S)
    165177
    V = conj.dot(U / S[None, :])
    166178
    else:
    @@ -171,7 +183,7 @@ def partial_svd(matrix, n_eigenvecs=None):
    171183
    # use dense form when sparse form will fail
    172184
    S, V = scipy.linalg.eigh(xTx.toarray())
    173185
    else:
    174-
    S, V = scipy.sparse.linalg.eigsh(xTx, k=n_eigenvecs, which='LM')
    186+
    S, V = scipy.sparse.linalg.eigsh(xTx, k=n_eigenvecs, which='LM', v0=v0)
    175187
    S = np.sqrt(S)
    176188
    U = matrix.dot(V / S[None, :])
    177189

    ‎tensorly/decomposition/_tucker.py

    Lines changed: 3 additions & 3 deletions
    Original file line numberDiff line numberDiff line change
    @@ -42,7 +42,7 @@ def partial_tucker(tensor, modes, rank=None, n_iter_max=100, init='svd', tol=10e
    4242
    4343
    Returns
    4444
    -------
    45-
    core : ndarray
    45+
    core : ndarray
    4646
    core tensor of the Tucker decomposition
    4747
    factors : ndarray list
    4848
    list of factors of the Tucker decomposition.
    @@ -74,7 +74,7 @@ def partial_tucker(tensor, modes, rank=None, n_iter_max=100, init='svd', tol=10e
    7474
    if init == 'svd':
    7575
    factors = []
    7676
    for index, mode in enumerate(modes):
    77-
    eigenvecs, _, _ = svd_fun(unfold(tensor, mode), n_eigenvecs=rank[index])
    77+
    eigenvecs, _, _ = svd_fun(unfold(tensor, mode), n_eigenvecs=rank[index], random_state=random_state)
    7878
    factors.append(eigenvecs)
    7979
    else:
    8080
    rng = check_random_state(random_state)
    @@ -87,7 +87,7 @@ def partial_tucker(tensor, modes, rank=None, n_iter_max=100, init='svd', tol=10e
    8787
    for iteration in range(n_iter_max):
    8888
    for index, mode in enumerate(modes):
    8989
    core_approximation = multi_mode_dot(tensor, factors, modes=modes, skip=index, transpose=True)
    90-
    eigenvecs, _, _ = svd_fun(unfold(core_approximation, mode), n_eigenvecs=rank[index])
    90+
    eigenvecs, _, _ = svd_fun(unfold(core_approximation, mode), n_eigenvecs=rank[index], random_state=random_state)
    9191
    factors[index] = eigenvecs
    9292

    9393
    core = multi_mode_dot(tensor, factors, modes=modes, transpose=True)

    ‎tensorly/decomposition/tests/test_tucker.py

    Lines changed: 8 additions & 1 deletion
    Original file line numberDiff line numberDiff line change
    @@ -3,7 +3,7 @@
    33
    from ...tucker_tensor import tucker_to_tensor
    44
    from ...tenalg import multi_mode_dot
    55
    from ...random import check_random_state
    6-
    from ...testing import assert_equal, assert_
    6+
    from ...testing import assert_equal, assert_, assert_array_equal
    77

    88

    99
    def test_partial_tucker():
    @@ -32,6 +32,13 @@ def test_partial_tucker():
    3232
    assert_equal(core.shape, [tensor.shape[0]]+ranks, err_msg="Core.shape={}, "
    3333
    "expected {}".format(core.shape, [tensor.shape[0]]+ranks))
    3434

    35+
    # Test random_state fixes the core and the factor matrices
    36+
    core1, factors1 = partial_tucker(tensor, modes=modes, rank=ranks, random_state=0)
    37+
    core2, factors2 = partial_tucker(tensor, modes=modes, rank=ranks, random_state=0)
    38+
    assert_array_equal(core1, core2)
    39+
    for factor1, factor2 in zip(factors1, factors2):
    40+
    assert_array_equal(factor1, factor2)
    41+
    3542

    3643
    def test_tucker():
    3744
    """Test for the Tucker decomposition"""

    ‎tensorly/tests/test_backend.py

    Lines changed: 9 additions & 0 deletions
    Original file line numberDiff line numberDiff line change
    @@ -411,6 +411,15 @@ def test_svd():
    411411
    true_rec_error = tl.sum((matrix - tl.dot(U, tl.reshape(S, (-1, 1))*V))**2)
    412412
    assert_(true_rec_error <= tol)
    413413

    414+
    # Test if partial_svd returns the same result for the same setting
    415+
    matrix = T.tensor(np.random.random((20, 5)))
    416+
    random_state = np.random.RandomState(0)
    417+
    U1, S1, V1 = tl.partial_svd(matrix, n_eigenvecs=2, random_state=random_state)
    418+
    U2, S2, V2 = tl.partial_svd(matrix, n_eigenvecs=2, random_state=0)
    419+
    assert_array_equal(U1, U2)
    420+
    assert_array_equal(S1, S2)
    421+
    assert_array_equal(V1, V2)
    422+
    414423

    415424
    def test_shape():
    416425
    A = T.arange(3*4*5)

    0 commit comments

    Comments
     (0)
    0