Skip to content

[MRG+1] Make csr row norms support fused types #6785

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 13 additions & 8 deletions sklearn/utils/sparsefuncs_fast.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -23,24 +23,29 @@ ctypedef np.float64_t DOUBLE

def csr_row_norms(X):
"""L2 norm of each row in CSR matrix X."""
if X.dtype != np.float32:
X = X.astype(np.float64)
return _csr_row_norms(X.data, X.shape, X.indices, X.indptr)


def _csr_row_norms(np.ndarray[floating, ndim=1, mode="c"] X_data,
shape,
np.ndarray[int, ndim=1, mode="c"] X_indices,
np.ndarray[int, ndim=1, mode="c"] X_indptr):
cdef:
unsigned int n_samples = X.shape[0]
unsigned int n_features = X.shape[1]
unsigned int n_samples = shape[0]
unsigned int n_features = shape[1]
np.ndarray[DOUBLE, ndim=1, mode="c"] norms
np.ndarray[DOUBLE, ndim=1, mode="c"] data
np.ndarray[int, ndim=1, mode="c"] indices = X.indices
np.ndarray[int, ndim=1, mode="c"] indptr = X.indptr

np.npy_intp i, j
double sum_

norms = np.zeros(n_samples, dtype=np.float64)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You might want to specify the type of norms here, i.e np.ndarray. There is some non-negligible overhead by not doing that. Btw, you can check that by compiling using the cython -a flag and checking the html.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it is just above, isn't it?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

hmm yeah.

data = np.asarray(X.data, dtype=np.float64) # might copy!

for i in range(n_samples):
sum_ = 0.0
for j in range(indptr[i], indptr[i + 1]):
sum_ += data[j] * data[j]
for j in range(X_indptr[i], X_indptr[i + 1]):
sum_ += X_data[j] * X_data[j]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you cast each entry to double as you multiply it, does it take a big runtime hit?

Copy link
Contributor Author

@yenchenlin yenchenlin May 19, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

About this, I've timed the csr_row_norms function with number of executions = 10000.

Original: 9.901501894
This PR: 9.89375782013
This PR: 13.0824899673

It seems that casting here won't become a bottleneck.

You can find my test script here.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably won't be.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My suggestion is that you can use float64 sum_ type and cast here to float64 so that you retain FP precision... if that's a real concern.

norms[i] = sum_

return norms
Expand Down
25 changes: 17 additions & 8 deletions sklearn/utils/tests/test_extmath.py
Original file line number Diff line number Diff line change
Expand Up @@ -148,14 +148,23 @@ def test_norm_squared_norm():

def test_row_norms():
X = np.random.RandomState(42).randn(100, 100)
sq_norm = (X ** 2).sum(axis=1)

assert_array_almost_equal(sq_norm, row_norms(X, squared=True), 5)
assert_array_almost_equal(np.sqrt(sq_norm), row_norms(X))

Xcsr = sparse.csr_matrix(X, dtype=np.float32)
assert_array_almost_equal(sq_norm, row_norms(Xcsr, squared=True), 5)
assert_array_almost_equal(np.sqrt(sq_norm), row_norms(Xcsr))
for dtype in (np.float32, np.float64):
if dtype is np.float32:
precision = 4
else:
precision = 5

X = X.astype(dtype)
sq_norm = (X ** 2).sum(axis=1)

assert_array_almost_equal(sq_norm, row_norms(X, squared=True),
precision)
assert_array_almost_equal(np.sqrt(sq_norm), row_norms(X), precision)

Xcsr = sparse.csr_matrix(X, dtype=dtype)
assert_array_almost_equal(sq_norm, row_norms(Xcsr, squared=True),
precision)
assert_array_almost_equal(np.sqrt(sq_norm), row_norms(Xcsr), precision)


def test_randomized_svd_low_rank_with_noise():
Expand Down