Consider this code:

import cr.sparse.data
import cr.sparse.lop
import jax.random
import jax.numpy as jnp
import matplotlib.pyplot as plt
import numpy as np

def to_array(op):
    return np.asarray(cr.sparse.lop.to_matrix(op))

def rand_orth_ref(key, n):
    Q, R = jnp.linalg.qr(jax.random.normal(key, (n, n)))
    return Q * jnp.sign(jnp.diag(R))  # this multiplication is important

def rand_orth_nonuniform(key, n):
    Q, R = jnp.linalg.qr(jax.random.normal(key, (n, n)))
    return Q

def rand_orth_nonuniform2(key, n):
    Q = rand_orth_ref(key, n)
    if Q[0, 0] > 0:
        Q = -Q  # Pretend the QR implementation returned -Q, -R instead (which it is allowed to do)
    return Q

def rand_orth_rows_dict(key, n):
    return to_array(cr.sparse.lop.random_orthonormal_rows_dict(key, n, n))

def rand_orth_onb_dict(key, n):
    return to_array(cr.sparse.lop.random_onb_dict(key, n))

def rand_orth_subspaces(key, n):
    Q, = cr.sparse.data.random_subspaces(key, n, n, 1)
    return np.asarray(Q)

def uniform_random(n, rand_orth):
    Qs = [rand_orth(key, 2) for key in jax.random.split(jax.random.PRNGKey(0), n)]
    xs, ys = np.array([Q[0] for Q in Qs]).T
    return np.arctan2(ys, xs) / (2 * np.pi) + .5

# ANGLES = uniform_random(10_000, rand_orth_ref)
# ANGLES = uniform_random(10_000, rand_orth_nonuniform)
# ANGLES = uniform_random(10_000, rand_orth_nonuniform2)
# ANGLES = uniform_random(10_000, rand_orth_rows_dict)
# ANGLES = uniform_random(10_000, rand_orth_onb_dict)
ANGLES = uniform_random(10_000, rand_orth_subspaces)

assert ANGLES.min() >= 0 and ANGLES.max() <= 1
print(ANGLES[ANGLES < .5].max(), ANGLES[ANGLES > .5].min())  # both should be close to 0.5
plt.hist(ANGLES, bins=20)
plt.xlim(0, 1);

All of the ways to generate ANGLES except for rand_orth_ref() are obviously non-uniform at least on my machine (the exact implementation of the QR decomposition generally depends on the system, but I don't know about JAX specifically). The print statement produces a result close to [0.25, 0.75] for me, i.e. in 2d half of the angles are never produced for the first basis vector.

The root cause is e.g. this code, which essentially implements sampling using the method rand_orth_nonuniform() above. The code is clearly intended to produce uniformly random samples, but it is missing the correction for signs as in rand_orth_ref(), which leaves the distribution dependent on how the QR decomposition happens to choose the arbitrary signs of the diagonal of R. Basically, QR = Q'R' where Q' = QD, R' = DR and D is any diagonal matrix with ±1 entries on the diagonal, and you can produce a non-uniform distribution by choosing D suitably, as in the rand_orth_nonuniform2() example. The choice of D is implementation-defined and hence arbitrary. The sign correction is mentioned e.g. in this Wikipedia article, which cites a Notices of the AMS article that goes into detail on the issue.

There are other places in the code using QR decomposition in the same way, at least here and here.

Hi @shailesh1729

I found the readme is somewhat cluttered.

Here are a few of my suggestions;

• Quickstart section is not necessary for this readme. I recommend removing it

• Please include the image about the sparse representation problem(at draft summary) below lines[14-25].

• Lines[27-31] can be replaced by For detailed documentation and usage, please visit https://cr-sparse.readthedocs.io/en/latest.

• Lines [34-51] looks redundant. I suggest removing it since you have mentioned the algorithms already in the introduction part and are present in the Examples section also. For more examples please direct the users to the example gallery https://cr-sparse.readthedocs.io/en/latest/gallery/index.html

• Line 57 is redundant. Remove line 57 and move line 56 to line 73

Thank you and let me know what you think about the suggestions,

This work provides implementations of algorithms in the paper: ALTERNATING DIRECTION ALGORITHMS FOR ℓ1-PROBLEMS IN COMPRESSIVE SENSING by Junfeng Yang and Yin Zhang.

Algorithms:

• Solve basis pursuit (1.3)
• Solve basis pursuit with inequality constraints (1.4)
• Solve basis pursuit denoising (1.5)
• Solve l_1/l_1 problem (1.6)

We also need to implement corresponding non-negative counterparts.

Evaluation:

• Recovery results from data corrupted by impulsive noise only
• Recovery results from data corrupted by white noise only
• Recovery results from data corrupted by both white and impulsive noise
• Relative error vs optimality for noiseless and noisy data

I feel that the following reference should be more visible in the online documentation:

In addition, please use the full reference for the paper, i.e., add author names, journal, year, vol, etc.

Hi @shailesh1729,

I have some suggestions about the paper. Let me know if you disagree with any of them:

Paper

Summary:

• The paragraphs between lines 14-61 can be under a separate heading, e.g.: 'Package Overview'.
• Brief introduction to problems of sparse representation / compressive sensing methods and the need for sparse recovery algorithms and their application domains for a broad and non-specialist audience. If you think it is out of the scope of this paper, refer the reader to the references where they can read the explanations.

Statement of need:

• This section could be used to give a general idea of what to expect from the package in terms of the concepts/problems mentioned in the Summary section, which I think are convincingly defined between lines 63 - 76. However, lines 77-83 are subjective. For example.

  • Line 77: There are some libraries (please cite them)
  • Line 81: There are several attempts (please cite them)
  • Line 83: Implementation appears to be somewhat complex (please explain).

• The rest of this section deals with JAX, which I think is beyond the scope of this project.

• The package has a broader scope, I suggest adding or naming some target groups.

Documentation/Readme

• Please give a brief introduction to sparse signal processing and the need for algorithms and the cr-sparse package. You can just copy this from the 'summary' section of the paper.

• Readme lines 17-40
  For ease of use, it would be very useful to classify these methods either by algorithms (sparse recovery/solvers) or by their applications (compressive sensing/sparse signal processing). Some examples of sparse recovery algorithm classification can be found here and here. IIf you feel that it is not feasible to provide the overview of the algorithms/methods classification in the readme file, I still recommend including it in the documentation and attaching the link in the readme file.

General

• Please provide some limitations of the package and future work. This would let the users know what you currently working on and what they cannot expect from the current version of the package.

Hello!

I hope you are doing well!

We are a security research team. Our tool automatically detected a vulnerability in this repository. We want to disclose it responsibly. GitHub has a feature called Private vulnerability reporting, which enables security research to privately disclose a vulnerability. Unfortunately, it is not enabled for this repository.

Can you enable it, so that we can report it?

Thanks in advance!

PS: you can read about how to enable private vulnerability reporting here: https://docs.github.com/en/code-security/security-advisories/repository-security-advisories/configuring-private-vulnerability-reporting-for-a-repository

I installed cr-sparse according to the installation instructions and tried following one of the tutorials, but ran into an Attribute error when trying to import cr-sparse. See below for the stack trace. I suspect there has been some update in Jax that isn't reflected yet in cr-sparse.

Package versions:
Jax 0.4.18
cr-sparse 0.3.2

System:
miniconda + VS Code
macOS Ventura 13.4, Apple M2

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
[/Users/joey/USC/lidar_group/shadowOQS/dynamical_shadow_tomography/cs_testing.ipynb](https://file+.vscode-resource.vscode-cdn.net/Users/joey/USC/lidar_group/shadowOQS/dynamical_shadow_tomography/cs_testing.ipynb) Cell 1 line 9

import matplotlib.pyplot as plt
# import cr.sparse as crs
--> import cr.sparse.dict as crdict
import cr.sparse.pursuit as pursuit

File [~/miniconda3/envs/sparse/lib/python3.11/site-packages/cr/sparse/__init__.py:37](https://file+.vscode-resource.vscode-cdn.net/Users/joey/USC/lidar_group/shadowOQS/dynamical_shadow_tomography/~/miniconda3/envs/sparse/lib/python3.11/site-packages/cr/sparse/__init__.py:37)
     31 # Evaluation Tools
     33 from cr.sparse._src.tools.performance import (
     34     RecoveryPerformance
     35 )
---> 37 from cr.sparse._src.tools.trials_at_fixed_m_n import (
     38     RecoveryTrialsAtFixed_M_N
     39 )

File [~/miniconda3/envs/sparse/lib/python3.11/site-packages/cr/sparse/_src/tools/trials_at_fixed_m_n.py:25](https://file+.vscode-resource.vscode-cdn.net/Users/joey/USC/lidar_group/shadowOQS/dynamical_shadow_tomography/~/miniconda3/envs/sparse/lib/python3.11/site-packages/cr/sparse/_src/tools/trials_at_fixed_m_n.py:25)
     22 import jax.numpy as jnp
     24 import cr.sparse as crs
---> 25 from cr.sparse import pursuit
     26 import cr.sparse.data as crdata
     27 import cr.sparse.dict as crdict

File [~/miniconda3/envs/sparse/lib/python3.11/site-packages/cr/sparse/pursuit/__init__.py:26](https://file+.vscode-resource.vscode-cdn.net/Users/joey/USC/lidar_group/shadowOQS/dynamical_shadow_tomography/~/miniconda3/envs/sparse/lib/python3.11/site-packages/cr/sparse/pursuit/__init__.py:26)
     18 # pylint: disable=W0611
     20 from cr.sparse._src.pursuit.util import (
     21     abs_max_idx,
     22     gram_chol_update,
     23     largest_indices
     24 )
---> 26 from cr.sparse._src.pursuit.defs import (
     27     RecoverySolution
     28 )

File [~/miniconda3/envs/sparse/lib/python3.11/site-packages/cr/sparse/_src/pursuit/defs.py:24](https://file+.vscode-resource.vscode-cdn.net/Users/joey/USC/lidar_group/shadowOQS/dynamical_shadow_tomography/~/miniconda3/envs/sparse/lib/python3.11/site-packages/cr/sparse/_src/pursuit/defs.py:24)
     20 from cr.nimble.dsp import build_signal_from_indices_and_values
     21 norm = jnp.linalg.norm
     23 @dataclass
---> 24 class SingleRecoverySolution:
     25     signals: jnp.DeviceArray = None
     26     representations : jnp.DeviceArray = None

File [~/miniconda3/envs/sparse/lib/python3.11/site-packages/cr/sparse/_src/pursuit/defs.py:25](https://file+.vscode-resource.vscode-cdn.net/Users/joey/USC/lidar_group/shadowOQS/dynamical_shadow_tomography/~/miniconda3/envs/sparse/lib/python3.11/site-packages/cr/sparse/_src/pursuit/defs.py:25), in SingleRecoverySolution()
     23 @dataclass
     24 class SingleRecoverySolution:
---> 25     signals: jnp.DeviceArray = None
     26     representations : jnp.DeviceArray = None
     27     residuals : jnp.DeviceArray =  None

File [~/miniconda3/envs/sparse/lib/python3.11/site-packages/jax/_src/deprecations.py:53](https://file+.vscode-resource.vscode-cdn.net/Users/joey/USC/lidar_group/shadowOQS/dynamical_shadow_tomography/~/miniconda3/envs/sparse/lib/python3.11/site-packages/jax/_src/deprecations.py:53), in deprecation_getattr.<locals>.getattr(name)
     51   warnings.warn(message, DeprecationWarning, stacklevel=2)
     52   return fn
---> 53 raise AttributeError(f"module {module!r} has no attribute {name!r}")

AttributeError: module 'jax.numpy' has no attribute 'DeviceArray'

• I wouldn't use attributive adjective such as "decent", as you used, in line 75 to describe the implementation in aaren/wavelets, or "old" in lines 76, 83.

• I think subjective claims, e.g., line 131 "we believe that our implementation is cleaner and simpler...", should not be included in the paper.