My environment is windows10 and install all package the project need.
Cython 0.29.20
numpy 1.19.0
scipy 1.5.1
scikit-learn 0.23.1
nose 1.3.7

I have been setup ok already ,and after successfully installing the ristretto library.
The unit tests can be run by: python setup.py test
When i run this code and show an error below:

(rnmf) C:\Users\tsai\Desktop\ristretto-master>python setup.py test
running test
WARNING: Testing via this command is deprecated and will be removed in a future version. Users looking for a generic test entry point independent of test runner are encouraged to use tox.
running egg_info
writing ristretto.egg-info\PKG-INFO
writing dependency_links to ristretto.egg-info\dependency_links.txt
writing requirements to ristretto.egg-info\requires.txt
writing top-level names to ristretto.egg-info\top_level.txt
reading manifest file 'ristretto.egg-info\SOURCES.txt'
reading manifest template 'MANIFEST.in'
warning: no files found matching '.c' under directory 'ristretto'
warning: no files found matching '.pyx' under directory 'ristretto'
warning: no files found matching 'README.md'
writing manifest file 'ristretto.egg-info\SOURCES.txt'
running build_ext
C:\Users\tsai\anaconda3\envs\rnmf\lib\site-packages\sklearn\utils\deprecation.py:143: FutureWarning: The sklearn.decomposition.cdnmf_fast module is deprecated in version 0.22 and will be removed in version 0.24. The corresponding classes / functions should instead be imported from sklearn.decomposition. Anything that cannot be imported from sklearn.decomposition is now part of the private API.
warnings.warn(message, FutureWarning)
C:\Users\tsai\anaconda3\envs\rnmf\lib\site-packages\sklearn\utils\deprecation.py:143: FutureWarning: The sklearn.decomposition.nmf module is deprecated in version 0.22 and will be removed in version 0.24. The corresponding classes / functions should instead be imported from sklearn.decomposition. Anything that cannot be imported from sklearn.decomposition is now part of the private API.
warnings.warn(message, FutureWarning)
test_sketches.test_random_axis_sample ... ok
test_sketches.test_random_gaussian_map ... ok
test_sketches.test_sparse_random_map ... ok
test_transforms.test_randomized_uniform_sampling ... ok
test_transforms.test_johnson_linderstrauss ... ok
test_transforms.test_sparse_johnson_linderstrauss ... ok
test_transforms.test_fast_johnson_linderstrauss ... ok
test_utils.test_orthonormalize ... ok
test_utils.test_perform_subspace_iterations ... ok
ristretto.tests.test_cur.test_compute_cur ... C:\Users\tsai\Desktop\ristretto-master\ristretto\utils.py:12: DeprecationWarning: Converting np.complex to a dtype is deprecated. The current result is complex128 which is not strictly correct.
if A.dtype == np.complexfloating:
ok
ristretto.tests.test_cur.test_compute_rcur ... ok
ristretto.tests.test_dmd.test_dmd ... ok
ristretto.tests.test_dmd.test_compute_rdmd ... ok
ristretto.tests.test_dmd.test_DMD ... ok
ristretto.tests.test_dmd.test_RDMD ... ok
ristretto.tests.test_eigen.test_compute_reigh_float64 ... ok
ristretto.tests.test_eigen.test_compute_reigh_complex128 ... ok
ristretto.tests.test_eigen.test_reig_nystroem_float64 ... ok
ristretto.tests.test_eigen.test_reig_nystroem_complex128 ... ok
ristretto.tests.test_eigen.test_reig_nystroem_col_float64 ... ok
ristretto.tests.test_eigen.test_reig_nystroem_col_complex128 ... ok
ristretto.tests.test_interp_decomp.test_id_col ... ok
ristretto.tests.test_interp_decomp.test_id_row ... ok
ristretto.tests.test_interp_decomp.test_rid_col ... ok
ristretto.tests.test_interp_decomp.test_rid_row ... ok
ristretto.tests.test_lu.test_compute_rlu_float64 ... ok
ristretto.tests.test_lu.test_compute_rlu_complex128 ... ok
ristretto.tests.test_nmf.test_nmf_fhals ... ERROR
ristretto.tests.test_nmf.test_rnmf_fhals ... ERROR
ristretto.tests.test_pca.test_spca ... ok
ristretto.tests.test_pca.test_robspca ... ok
ristretto.tests.test_pca.test_rspca ... ok
ristretto.tests.test_qb.test_rqb_float64 ... ok
ristretto.tests.test_qb.test_rqb_complex128 ... ok
ristretto.tests.test_qb.test_rqb_block_float64 ... ok
ristretto.tests.test_qb.test_rqb_block_wide_float64 ... ok
ristretto.tests.test_qb.test_rqb_block_complex128 ... ok
ristretto.tests.test_qb.test_rqb_block_wide_complex128 ... ok
ristretto.tests.test_svd.test_compute_rsvd_float64 ... ok
ristretto.tests.test_svd.test_compute_rsvd_complex128 ... ok

======================================================================
ERROR: ristretto.tests.test_nmf.test_nmf_fhals

Traceback (most recent call last):
File "C:\Users\tsai\anaconda3\envs\rnmf\lib\site-packages\nose\case.py", line 198, in runTest
self.test(*self.arg)
File "C:\Users\tsai\Desktop\ristretto-master\ristretto\tests\test_nmf.py", line 15, in test_nmf_fhals
W, H = compute_nmf(Anoisy, rank=10)
File "C:\Users\tsai\Desktop\ristretto-master\ristretto\nmf.py", line 154, in compute_nmf
violation = _update_cdnmf_fast(Ht, WtW, AtW, permutation)
File "sklearn\decomposition_cdnmf_fast.pyx", line 13, in sklearn.decomposition._cdnmf_fast._update_cdnmf_fast
ValueError: Buffer dtype mismatch, expected 'Py_ssize_t' but got 'long'

======================================================================
ERROR: ristretto.tests.test_nmf.test_rnmf_fhals

Traceback (most recent call last):
File "C:\Users\tsai\anaconda3\envs\rnmf\lib\site-packages\nose\case.py", line 198, in runTest
self.test(*self.arg)
File "C:\Users\tsai\Desktop\ristretto-master\ristretto\tests\test_nmf.py", line 26, in test_rnmf_fhals
W, H = compute_rnmf(Anoisy, rank=10)
File "C:\Users\tsai\Desktop\ristretto-master\ristretto\nmf.py", line 335, in compute_rnmf
violation = _update_cdnmf_fast(Ht, WtW, BtW, permutation)
File "sklearn\decomposition_cdnmf_fast.pyx", line 13, in sklearn.decomposition._cdnmf_fast._update_cdnmf_fast
ValueError: Buffer dtype mismatch, expected 'Py_ssize_t' but got 'long'

Ran 40 tests in 0.850s

FAILED (errors=2)
Test failed: <unittest.runner.TextTestResult run=40 errors=2 failures=0>
error: Test failed: <unittest.runner.TextTestResult run=40 errors=2 failures=0>

Currently the tests all live in the same file tests/tests.py. I think it would be better if we split up the testing into directories corresponding to the modules of the package (ie dmd, mf, ...)

Move to a more object oriented framework for easier access of attributes of our algorithms.

• add classes compatible with scikit-learn estimators
• add compute_ prefix to current functions

This should end up looking something like:

ristretto
|
---> dmd
        |
         ----> compute_dmd
         ----> DMD

This looks like a great project, but I've had a couple problems trying to install and use ristretto based on the documentation. I'm on macOS 10.13.4. First, I tried to simply pip install it, but it failed because of a Cython dependency issue. I installed Cython, and the installation still failed, throwing errors like,

gcc -Wno-unused-result -Wsign-compare -Wunreachable-code -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -I/Users/gwg/miniconda3/envs/ele538b/include -arch x86_64 -I/Users/gwg/miniconda3/envs/ele538b/include -arch x86_64 -I/Users/gwg/miniconda3/envs/ele538b/include/python3.7m -c ristretto/nmf/_fhals_update.c -o build/temp.macosx-10.7-x86_64-3.7/ristretto/nmf/_fhals_update.o
ristretto/nmf/_fhals_update.c:16255:21: error: no member named 'exc_type' in 'struct _ts'
    *type = tstate->exc_type;
            ~~~~~~  ^

Eventually I got a fatal error: too many errors emitted, and the installation process crashed.

So then I cloned the directory and ran setup.py instead, which almost worked, but I also had to install sckit-learn before I could import it. Now I can import the module, but I get errors when importing specific functions, e.g.:

ImportError: cannot import name 'rlu' from 'ristretto.lu' (/Users/gwg/miniconda3/envs/ele538b/lib/python3.7/site-packages/ristretto-0.1.2-py3.7.egg/ristretto/lu.py)

I apologize for the broad request, but could the maintainers please clean up the installation process or advise if I'm doing something wrong? Thanks.

Replace print() functions within algorithms with either:

• logging.<log_level>()
• warnings.<warn_level>()
  where appropriate, specified by user verbosity or verbose keyword arguments.

I want to propose a couple changes to the Ristretto API. I'm happy to try make these changes and create a pull request if the main contributors think this is a good idea.

1. Rather than name methods compute_rx where x is the name of the non-randomized method, what about just rx? For example, rsvd instead of compute_rsvd, rlu instead of compute_rlu. This would match libraries like NumPy, but with an r prefix to indicate it is a randomized version.

2. Rather than importing from modules named ristretto.x where x is the name of the non-randomized method, import directly from ristretto. For example, rather than importing from ristretto.cur, import directly from ristretto.

As an example of combining these two changes, instead of,

>>> from ristretto.svd import compute_rsvd

we would have,

>>> from ristretto import rsvd

I cannot import any modules or functions from the package:

from ristretto.svd import rsvd
Traceback (most recent call last):

  File "<ipython-input-2-c45526633c15>", line 1, in <module>
    from ristretto.svd import rsvd

ModuleNotFoundError: No module named 'ristretto.svd'

Also, I find it more handy, if we can do:

import ristretto as ro

so that I can call any function via

ro.rsvd(), etc.

Best,
B

Add support for arrays of type scipy.sparse.[coo_matrix, csr_matrix, ...].

This has a lot of potential gotchya's, however it would be nice to add support
for them in the algorithms in which it is possible to do so.

I test the accuracy of ristretto.pca.compute_spca by setting:

• alpha = 0. No L1 regularization.
• n_components = number of features. No dimension reduction.
• tol = 1e-32, max_iter=1e6. A large number of iterations.

I assume the result should be same as a standard PCA over matrix X.

Definition for relative error is: Error[i, j] = 100 * Abs(1 - X_{reconstructed}[i, j] / X_actual[i, j]).
X has no zeros or very small numbers.

However, reconstruction error can be large for some elements. I guess as the covariance matrix becomes more sparse, the spikes in the relative error becomes taller. But I am not sure.

import ristretto.pca
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import cholesky


def make_covariance_matrix(n_feat, zero_prob=0.99):
    shape = (n_feat, n_feat)
    cov_mat = np.random.uniform(0, 1, shape)
    cov_mat *= np.random.choice(2, shape, p=[zero_prob, 1 - zero_prob])
    np.fill_diagonal(cov_mat, 1)
    cov_mat = cov_mat.dot(cov_mat.T)
    return cov_mat


def make_independent_normal_dists(n_feat, n_obs):
    return np.random.normal(0, 1, (n_feat, n_obs))


def make_correlated_normal_dists(n_obs, n_feat, cov_mat=None):
    """
    https://scipy-cookbook.readthedocs.io/items/CorrelatedRandomSamples.html
    """
    if cov_mat is None:
        cov_mat = make_covariance_matrix(n_feat)
    c = cholesky(cov_mat, lower=True)
    ind_norm = make_independent_normal_dists(n_feat, n_obs)
    return np.dot(c, ind_norm).T, cov_mat


def plot_covariance_mat(cov_mat):
    plt.title('covariance matrix')
    plt.imshow(cov_mat)
    plt.colorbar()
    plt.xlabel('feature index')
    plt.ylabel('feature index')
    plt.show()
    plt.close()


def plot_filtered_error(rel_error, percentile_thres=99.5):
    error_1d = np.ravel(rel_error)
    threshold = np.percentile(rel_error, percentile_thres)
    error_filtered = error_1d[error_1d<threshold]
    plt.title('Histogram of Relative Error below %f quantile' %
              percentile_thres)
    plt.hist(np.ravel(error_filtered), bins=100)
    plt.xlabel('Relative Error')
    plt.ylabel('Count')
    plt.show()
    plt.close()


def main():
    n_obs = 5000  # Number of observations
    n_feat = 50  # Number of features
    n_components = 50  # Number of principal components.

    # Parameters for standard PCA (not sparse).
    sparsity_alpha = 0
    robust = False
    ridge_beta = 0.1

    # Create test data with zero mean along each column
    X, cov_mat = make_correlated_normal_dists(n_obs // 2, n_feat)
    X -= np.mean(X, axis=0)  # Centering each column by its mean.
    X /= np.std(X, axis=0)
    X = np.vstack([X - 4, X + 4])
    # np.random.shuffle(X)  # shuffle rows of X
    plot_covariance_mat(cov_mat)

    # Sparse PCA
    B, A, eigen_values, obj = ristretto.pca.compute_spca(
        X, n_components=n_components, alpha=sparsity_alpha, robust=robust,
        beta=ridge_beta, tol=1e-32, max_iter=1e6)

    # Reconstruction of X by (XB)A^{T}
    recons_x = np.matmul(np.matmul(X, B), np.transpose(A))

    # Calculate relative error
    delta = X - recons_x
    rel_error = np.abs(delta / X) * 100

    # Plot histogram of relative error
    # filter relative error because some are crazy.
    plot_filtered_error(rel_error)

    # Plot relative error
    plt.title('Relative Reconstruction Error (capped at 200%)')
    plt.imshow(rel_error, vmin=0, vmax=200, aspect='auto')
    plt.colorbar()
    plt.xlabel('feature index')
    plt.ylabel('observation index')
    plt.show()
    plt.close()


main()

Adding a random_state keyword argument to the algorithms make it easier for
the user to reproduce results for testing or validation.

RSPCA returns a bug at fit time, see the below example:

from ristretto.pca import RSPCA
import numpy as np

X = np.random.RandomState(42).randn(20, 20)
rspca = RSPCA()
rspca.fit(X)

returns:

TypeError: compute_rspca() got an unexpected keyword argument 'gamma'

I am getting an error message of TypeError: dmd() got an unexpected keyword argument 'rank' when I run the code test_dmd.py.

Standardize the behavior of the algorithms when n > m. Choices include:

• 'flipped' (current): initially transpose input array A, then transpose
  the returned results at the end.
• 'axis': compute sketch function on the row space and returning A.dot(S)
  instead of S.T.dot(A) for QB decomposition.

Maybe, for performance we should change A.dot(B) to sci.dot(A, B). Are the following results reproducible?

import numpy as np
import scipy as sci
import timeit

A = sci.random.rand(1500,1500)
A = np.float32(A)

t0 = timeit.default_timer()
for i in range(100):
    C = A.dot(A.T)
print('Time C = A.dot(A): %s' %(timeit.default_timer()  - t0) )


t0 = timeit.default_timer()
for i in range(100):
    C = np.dot(A, A.T)
print('Time C = np.dot(A, A): %s' %(timeit.default_timer()  - t0) )


t0 = timeit.default_timer()
for i in range(100):
    C = sci.dot(A, A.T)
print('Time C = sci.dot(A, A): %s' %(timeit.default_timer()  - t0) )


t0 = timeit.default_timer()
C = sci.empty_like(C)
for i in range(100):
    sci.dot(A, A.T, out = C)
print('Time sci.dot(A, A, out = C): %s' %(timeit.default_timer()  - t0) )

Results are:

Time C = A.dot(A.T): 6.120499561999168
Time C = np.dot(A, A.T): 5.945316646997526
Time C = sci.dot(A, A.T): 5.941789857002732
Time sci.dot(A, A, out = C): 5.817408093993436

Best,
Ben

This blocked version can compute rQB of hdf5 datasets too large to fit in memory, and can be incorporated into algorithms which use the rQB (rSVD, rDMD, ...)

Currently the entire package is 'flattened' meaning that you can import any function
in any module of the package directly by:

from ristretto import <function>

This is more or less equivalent to writing the entire package in a single file,
and may be undesirable especially because of namespace confusion.

As of now, from ristretto import dmd imports the function dmd
from the dmd module within the package, not the dmd module. In addition
to being confusing, this prevents a user from importing the dmd module at all,
ie:

from ristretto import dmd
...
standard_result = dmd.dmd(A)
random_result = dmd.rdmd(A)

To accomplish the above, the user as of now must:

import ristretto as ro
...
standard_result = ro.dmd(A)
random_result = ro.rdmd(A)

Which is fine if the package only has a few modules containing a couple functions
each, but it becomes much less manageable as modules are added.

In the end this comes down to the following:

1. Will the package only contain a handful of matrix factorization functions
   or
2. Will modules be added to the package that add new or different sets of methods
   (i.e. randomized graphical methods like randomized min cut or randomized regression
   techniques). Additionally, will the current modules be extended to include different
   computational paradigms (i.e. parallel methods or limited precision methods).

If 1, the flattening is fine, but if any part of 2 is wanted in the future,
I think it is best that we spend some time rethinking the package import structure.

This way we can check all current functions, and future merges are covered with tests.

Standardize the set of sampling distributions. As of now we use the following:

• 'guassian' or 'normal': np.random.standard_normal
• 'uniform': np.random.uniform(-1, 1)
• 'None': np.random.rand
  And also from csvd
• 'spixel': random sampling of rows
• 'sparse': np.sparse.rand

After we agree on the set of sampling distributions, we can incorporate them into
the sketch function to reduce code repetition and limit the set of distributions
on an algorithm by algorithm basis if needed (i.e. if rsvd and rdmd
should not have the same set of distributions available).

All observations form a matrix X where each row is one observation and each column is one feature.
Each observation also has a weight w_i, which means the observation has occurred for w_i times.
For example, a weight of 2 means the same observation actually occurs for twice.

Can I do weighted sparse PCA by the following method?

1. Center each column of X by weighted average with w_i as weights.
2. Form matrix Y by multiplying each row of X by square root of w_i.
3. Calculate (robust) sparse PCA of Y with ristretto, which provides matrices A and B.
4. The matrix X is approximately XBA^T, where A^T is transpose of A.

Thank you.

For each combination of n_components and alpha, I run the following for 6 times.
After that I plot eigen_values as a scree plot.

B, A, eigen_values, obj = ristretto.pca.compute_rspca(
    X, n_components=n_components, n_subspace=3, alpha=alpha, robust=robust)

For a small alpha and robust=True, the scree plot looks normal.
The horizontal axis is the component index, the vertical axis is the median eigenvalue over the 6 runs.
The color of the line corresponds to n_components.

However, for robust=False or higher values of alpha, I get a different kind of scree plot:

While nearly all of the plots shows the 3rd and higher components all have an eigenvalue close to 0, I don't quite believe the system only has 2 components.

Why do I get the 2nd kind of scree plot?

The last PyPI release of Ristretto was version 0.1.2 in August of 2017. Since then there have been several changes to Ristretto that might warrant a new release. To help with moving a potential release along, I provide an ultra-short summary of release steps below. (The first five-or-so times I deployed to PyPI I went through a cycle of forgetting and rediscovering these steps.)

The bare basics of the new release would consist of three steps. First, update

Beyond the three steps there, you can draft a GitHub release to record release notes. See https://github.com/cvxgrp/cvxpy/releases/tag/v1.1.10 for an example GitHub release. When making an API-breaking release you would want to update the web documentation. However, since the web docs already suggest installing from GitHub, my guess is the web docs are already compatible with a would-be version 0.1.3 or 0.2.

If you're trying to decide between releasing under 0.1.3 or 0.2, then I'd recommend looking at the semantic versioning guidelines https://semver.org/. My guess is that releasing under 0.2 would be appropriate.

I have a little error about LinAlgError exception in eigen decompositions like compute_reigh_nystroem.
The error message is: NameError: name 'LinAlgError' is not defined

I think it should be

try: blabla except linalg.LinAlgError: blabla

instead of
try: blabla except LinAlgError: blabla

I'm working in a linux container using docker in Windows 10.

• Python 3.7.6
• numpy 1.18.5
• scipy 1.4.1
  Is it necessary give more info about my hardware/software?

OffTopic: Is it repository unmantained?