Skip to content
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

[FEA] Lanczos solver v2 #2481

Merged
merged 38 commits into from
Nov 20, 2024
Merged
Changes from 1 commit
Commits
Show all changes
38 commits
Select commit Hold shift + click to select a range
2c68742
init
aamijar Aug 17, 2024
0f5bdcb
benchmarking lanczos working
aamijar Aug 20, 2024
5ee60cb
format style
aamijar Aug 21, 2024
fa222c1
eigsh pylibraft api
aamijar Aug 21, 2024
7240937
gtest
aamijar Aug 21, 2024
a4cdc7a
clean up code
aamijar Aug 21, 2024
2be393a
update gtest rng seed
aamijar Aug 22, 2024
4f37b5c
update gtest edge case
aamijar Aug 22, 2024
79ce3f1
update gtest clean
aamijar Aug 22, 2024
ceb1d7a
resolving pr comments
aamijar Aug 23, 2024
14b7266
resolving pr comments
aamijar Aug 23, 2024
c564352
resolving pr comments
aamijar Aug 23, 2024
aba9c4e
resolving pr comments
aamijar Aug 24, 2024
ba16aca
resolving pr comments
aamijar Aug 24, 2024
74908f2
resolving pr comments
aamijar Aug 24, 2024
7b31108
resolving pr comments
aamijar Aug 27, 2024
2cdcc66
resolving pr comments
aamijar Aug 28, 2024
5aa1df2
Merge branch 'branch-24.10' into lanczos-solver-new
aamijar Aug 30, 2024
e4afa2d
Merge branch 'branch-24.10' into lanczos-solver-new
aamijar Aug 30, 2024
1160722
resolving pr comments
aamijar Sep 9, 2024
e473728
resolving pr comments
aamijar Sep 9, 2024
cc22a39
resolving pr comments
aamijar Sep 9, 2024
8767c6a
resolving pr comments
aamijar Sep 9, 2024
a3809eb
resolving pr comments
aamijar Sep 9, 2024
d4b4955
resolving pr comments
aamijar Sep 9, 2024
2f317b8
Add python test
lowener Oct 30, 2024
428ead4
Resolve PR comments
lowener Nov 1, 2024
985db38
Merge branch 'branch-24.12' into lanczos-solver-new
lowener Nov 1, 2024
2a14b26
Fix style
lowener Nov 1, 2024
9bc731c
Fix doc
lowener Nov 2, 2024
b870aac
Fix merge
lowener Nov 4, 2024
288ea20
Improve detail namespace usage, fix documentation
lowener Nov 5, 2024
80dafb9
Add python doc
lowener Nov 11, 2024
f8078e1
Consider v0 as optional in c++
lowener Nov 14, 2024
16459c7
Merge branch 'branch-24.12' into lanczos-solver-new
cjnolet Nov 16, 2024
45ee43f
Merge branch 'branch-24.12' into lanczos-solver-new
lowener Nov 18, 2024
4d0d4e8
Add issue number in code
lowener Nov 18, 2024
210af0a
Use new lanczos solver in raft::sparse::spectral
lowener Nov 19, 2024
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
Prev Previous commit
Next Next commit
Add python test
lowener committed Oct 30, 2024
commit 2f317b8c91093d00a8b723fe84df80cfb6df3827
9 changes: 4 additions & 5 deletions python/pylibraft/pylibraft/solver/lanczos.pyx
Original file line number Diff line number Diff line change
@@ -111,10 +111,9 @@ def eigsh(A, k=6, v0=None, ncv=None, maxiter=None,
with corresponding eigenvectors ``x``.

Args:
a (ndarray, spmatrix or LinearOperator): A symmetric square matrix with
dimension ``(n, n)``. ``a`` must :class:`cupy.ndarray`,
:class:`cupyx.scipy.sparse.spmatrix` or
:class:`cupyx.scipy.sparse.linalg.LinearOperator`.
a (spmatrix): A symmetric square sparse CSR matrix with
dimension ``(n, n)``. ``a`` must be of type
:class:`cupyx.scipy.sparse._csr.csr_matrix`
k (int): The number of eigenvalues and eigenvectors to compute. Must be
``1 <= k < n``.
v0 (ndarray): Starting vector for iteration. If ``None``, a random
@@ -175,7 +174,7 @@ def eigsh(A, k=6, v0=None, ncv=None, maxiter=None,
tol = np.finfo(ValueType).eps

if v0 is None:
rng = np.random.default_rng(seed)
rng = cp.random.default_rng(seed)
v0 = rng.random((N,)).astype(vals.dtype)

v0 = cai_wrapper(v0)
136 changes: 136 additions & 0 deletions python/pylibraft/pylibraft/test/test_solver.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,136 @@
# Copyright (c) 2024, NVIDIA CORPORATION.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#

import numpy
import pytest
from pylibraft.solver import eigsh
import cupy
from cupyx.scipy import sparse
import cupyx.scipy.sparse.linalg # NOQA

def shaped_random(
shape, xp=cupy, dtype=numpy.float32, scale=10, seed=0, order='C'):
"""Returns an array filled with random values.

Args:
shape(tuple): Shape of returned ndarray.
xp(numpy or cupy): Array module to use.
dtype(dtype): Dtype of returned ndarray.
scale(float): Scaling factor of elements.
seed(int): Random seed.

Returns:
numpy.ndarray or cupy.ndarray: The array with
given shape, array module,

If ``dtype`` is ``numpy.bool_``, the elements are
independently drawn from ``True`` and ``False``
with same probabilities.
Otherwise, the array is filled with samples
independently and identically drawn
from uniform distribution over :math:`[0, scale)`
with specified dtype.
"""
numpy.random.seed(seed)
dtype = numpy.dtype(dtype)
if dtype == '?':
a = numpy.random.randint(2, size=shape)
elif dtype.kind == 'c':
a = numpy.random.rand(*shape) + 1j * numpy.random.rand(*shape)
a *= scale
else:
a = numpy.random.rand(*shape) * scale
return xp.asarray(a, dtype=dtype, order=order)


class TestEigsh:
n = 30
density = 0.33
tol = {numpy.float32: 1e-5, numpy.complex64: 1e-5, 'default': 1e-12}
res_tol = {'f': 1e-5, 'd': 1e-12}
return_eigenvectors = True

def _make_matrix(self, dtype, xp):
shape = (self.n, self.n)
a = shaped_random(shape, xp, dtype=dtype)
mask = shaped_random(shape, xp, dtype='f', scale=1)
a[mask > self.density] = 0
a = a * a.conj().T
return a

def _test_eigsh(self, a, k, xp, sp):
expected_ret = sp.linalg.eigsh(a, k=k, return_eigenvectors=self.return_eigenvectors)
actual_ret = eigsh(a, k=k)
if self.return_eigenvectors:
w, x = actual_ret
exp_w, _ = expected_ret
# Check the residuals to see if eigenvectors are correct.
ax_xw = a @ x - xp.multiply(x, w.reshape(1, k))
res = xp.linalg.norm(ax_xw) / xp.linalg.norm(w)
tol = self.res_tol[numpy.dtype(a.dtype).char.lower()]
assert (res < tol)
else:
w = actual_ret
exp_w = expected_ret
w = xp.sort(w)
cupy.allclose(w, exp_w, rtol=tol, atol=tol)

@pytest.mark.parametrize('format', ['csr']) # , 'csc', 'coo'])
@pytest.mark.parametrize('k', [3, 6, 12])
@pytest.mark.parametrize('dtype', ['f', 'd'])
def test_sparse(self, format, k, dtype, xp=cupy, sp=sparse):
if format == 'csc':
pytest.xfail('may be buggy') # trans=True

a = self._make_matrix(dtype, xp)
a = sp.coo_matrix(a).asformat(format)
return self._test_eigsh(a, k, xp, sp)

def test_invalid(self):
xp, sp = cupy, sparse
a = xp.diag(xp.ones((self.n, ), dtype='f'))
with pytest.raises(ValueError):
sp.linalg.eigsh(xp.ones((2, 1), dtype='f'))
with pytest.raises(ValueError):
sp.linalg.eigsh(a, k=0)
a = xp.diag(xp.ones((self.n, ), dtype='f'))
with pytest.raises(ValueError):
sp.linalg.eigsh(xp.ones((1,), dtype='f'))
with pytest.raises(TypeError):
sp.linalg.eigsh(xp.ones((2, 2), dtype='i'))
with pytest.raises(ValueError):
sp.linalg.eigsh(a, k=self.n)

def test_starting_vector(self):
n = 100

# Make symmetric matrix
aux = self._make_matrix('f', cupy)
aux = sparse.coo_matrix(aux).asformat('csr')
matrix = (aux + aux.T) / 2.0

# Find reference eigenvector
ew, ev = eigsh(matrix, k=1)
v = ev[:, 0]

# Obtain non-converged eigenvector from random initial guess.
ew_aux, ev_aux = eigsh(matrix, k=1, ncv=1, maxiter=0)
v_aux = cupy.copysign(ev_aux[:, 0], v)

# Obtain eigenvector using known eigenvector as initial guess.
ew_v0, ev_v0 = eigsh(matrix, k=1, v0=v.copy(), ncv=1, maxiter=0)
v_v0 = cupy.copysign(ev_v0[:, 0], v)

assert cupy.linalg.norm(v - v_v0) < cupy.linalg.norm(v - v_aux)