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

add linalg.lstsq #2165

Merged
merged 9 commits into from May 17, 2019
@@ -21,6 +21,7 @@
from cupy.linalg.eigenvalue import eigvalsh # NOQA

from cupy.linalg.solve import inv # NOQA
from cupy.linalg.solve import lstsq # NOQA
from cupy.linalg.solve import pinv # NOQA
from cupy.linalg.solve import solve # NOQA
from cupy.linalg.solve import tensorinv # NOQA
@@ -166,7 +166,61 @@ def tensorsolve(a, b, axes=None):
return result.reshape(oldshape)


# TODO(okuta): Implement lstsq
def lstsq(a, b, rcond=1e-15):
"""Return the least-squares solution to a linear matrix equation.
Solves the equation `a x = b` by computing a vector `x` that
minimizes the Euclidean 2-norm `|| b - a x ||^2`. The equation may
be under-, well-, or over- determined (i.e., the number of
linearly independent rows of `a` can be less than, equal to, or
greater than its number of linearly independent columns). If `a`
is square and of full rank, then `x` (but for round-off error) is
the "exact" solution of the equation.
Args:
a (cupy.ndarray): "Coefficient" matrix with dimension ``(M, N)``
b (cupy.ndarray): "Dependent variable" values with dimension ``(M,)``
or ``(M, K)``
rcond (float): Cutoff parameter for small singular values.
For stability it computes the largest singular value denoted by
``s``, and sets all singular values smaller than ``s`` to zero.
Returns:
tuple:
A tuple of ``(x, residuals, rank, s)``. Note ``x`` is the
least-squares solution with shape ``(N,)`` or ``(N, K)`` depending
if ``b`` was two-dimensional. The sums of ``residuals`` is the
squared Euclidean 2-norm for each column in b - a*x. The
``residuals`` is an empty array if the rank of a is < N or M <= N,
but iff b is 1-dimensional, this is a (1,) shape array, Otherwise
the shape is (K,). The ``rank`` of matrix ``a`` is an integer. The
singular values of ``a`` are ``s``.
.. seealso:: :func:`numpy.linalg.lstsq`
"""
m, n = a.shape[-2:]
This conversation was marked as resolved by asi1024

This comment has been minimized.

Copy link
@asi1024

asi1024 May 7, 2019

Member

Could you check dimensions of the given array and raise linalg.LinAlgError?

This comment has been minimized.

Copy link
@cjekel

cjekel May 7, 2019

Author Contributor

I added checks for dimensions on a and b, as well as compatibility between a and b.

    util._assert_cupy_array(a, b)
    util._assert_rank2(a)
    if b.ndim > 2:
        raise linalg.LinAlgError('{}-dimensional array given. Array must be at'
                                 ' most two-dimensional'.format(b.ndim))
    m, n = a.shape[-2:]
    m2 = b.shape[0]
    if m != m2:
        raise linalg.LinAlgError('Incompatible dimensions')
u, s, vt = cupy.linalg.svd(a, full_matrices=False)
cutoff = rcond * s.max()
s1 = 1 / s
sing_vals = s <= cutoff
s1[sing_vals] = 0
rank = s.size - sing_vals.sum()
if b.ndim > 1:
s1 = cupy.repeat(s1.reshape(-1, 1), b.shape[1], axis=1)
z = core.dot(u.transpose(), b) * s1
x = core.dot(vt.transpose(), z)
if rank != n or m <= n:
resids = cupy.array([], dtype=a.dtype)
elif b.ndim > 1:
k = b.shape[1]
resids = cupy.zeros(k, dtype=a.dtype)
for i in range(k):
e = b[:, i] - core.dot(a, x[:, i])
resids[i] = core.dot(e.T, e)
This conversation was marked as resolved by asi1024

This comment has been minimized.

Copy link
@asi1024

asi1024 May 7, 2019

Member

This for-loop can be avoidable like as follows:

elif b.ndim > 1:
    e = b - core.dot(a, x)
    resids = cupy.sum(cupy.square(e), axis=0)

This comment has been minimized.

Copy link
@cjekel

cjekel May 7, 2019

Author Contributor

Thanks! I was wondering how to do this better.

else:
e = b - core.dot(a, x)
resids = core.dot(e.T, e).reshape(-1)
return x, resids, rank, s


def inv(a):
@@ -148,6 +148,65 @@ def test_invalid_shape(self):
self.check_shape((4, 3, 2, 1), rcond=0.1)


@unittest.skipUnless(
cuda.cusolver_enabled, 'Only cusolver in CUDA 8.0 is supported')
@testing.gpu
class TestLstsq(unittest.TestCase):

@testing.for_float_dtypes(no_float16=True)
@condition.retry(10)
def check_x(self, a_shape, b_shape, rcond, dtype):
a_cpu = numpy.random.randint(0, 10, size=a_shape).astype(dtype)
b_cpu = numpy.random.randint(0, 10, size=b_shape).astype(dtype)
a_gpu = cupy.asarray(a_cpu)
b_gpu = cupy.asarray(b_cpu)
a_gpu_copy = a_gpu.copy()
b_gpu_copy = b_gpu.copy()
result_cpu, resids_cpu, rank_cpu, s_cpu = numpy.linalg.lstsq(a_cpu,
b_cpu,
rcond=rcond) # noqa E501
result_gpu, resids_gpu, rank_gpu, s_gpu = cupy.linalg.lstsq(a_gpu,
b_gpu,
rcond=rcond) # noqa E501

self.assertEqual(result_cpu.dtype, result_gpu.dtype)
cupy.testing.assert_allclose(result_cpu, result_gpu, atol=1e-3)
cupy.testing.assert_allclose(resids_cpu, resids_gpu, atol=1e-3)
cupy.testing.assert_allclose(rank_cpu, rank_gpu, atol=1e-3)
cupy.testing.assert_allclose(s_cpu, s_gpu, atol=1e-3)
cupy.testing.assert_array_equal(a_gpu_copy, a_gpu)
cupy.testing.assert_array_equal(b_gpu_copy, b_gpu)

This comment has been minimized.

Copy link
@asi1024

asi1024 May 7, 2019

Member

Does not this implementation ensure that the result of cupy.lstsq is equal to that of numpy.lstsq? In other words, it ensures only ||ax-b||^2 are equal to each other?
If so, could you fix to check if ||ax-b||^2 is close?

This comment has been minimized.

Copy link
@cjekel

cjekel May 7, 2019

Author Contributor

I check if ||ax-b||^2 is close with cupy.testing.assert_allclose(resids_cpu, resids_gpu, atol=1e-3)

I'm not sure about atol=1e-3, I picked this because it was the same value for tests of cupy.linalg.pinv which also uses svd.

This comment has been minimized.

Copy link
@asi1024

asi1024 May 8, 2019

Member

Does not this implementation ensure that the result (the first element of the return values) of cupy.lstsq is equal to that of numpy.lstsq?

This comment has been minimized.

Copy link
@cjekel

cjekel May 8, 2019

Author Contributor

The implementation ensures that each element in the result is close, but not equal. They won't be equal in a programming, because the backends are different.

Edit: To clarify, the implementation checks that all elements in the results of numpy are close to cupy.

This comment has been minimized.

Copy link
@asi1024

asi1024 May 8, 2019

Member

Why @condition.retry(10) at L157 is needed? One of the result value is often distant from the expected value?

This comment has been minimized.

Copy link
@cjekel

cjekel May 8, 2019

Author Contributor

Why @condition.retry(10) at L157 is needed? One of the result value is often distant from the expected value?

Long response #2165 (comment)


def check_shape(self, a_shape, b_shape):
This conversation was marked as resolved by asi1024

This comment has been minimized.

Copy link
@asi1024

asi1024 May 7, 2019

Member

Why not clarifying that this method takes invalid shapes?

Suggested change
def check_shape(self, a_shape, b_shape):
def check_invalid_shape(self, a_shape, b_shape):

This comment has been minimized.

Copy link
@cjekel

cjekel May 7, 2019

Author Contributor

I also changed def check_x to def check_lstsq_solution.

a = cupy.random.rand(*a_shape)
b = cupy.random.rand(*b_shape)
with self.assertRaises(ValueError):
cupy.linalg.lstsq(a, b)

def test_lstsq(self):
# Test skinny
self.check_x((10, 3), (10, ), rcond=1e-15)
self.check_x((10, 3), (10, 2), rcond=1e-15)
# Test square
self.check_x((4, 4), (4, ), rcond=1e-15)
self.check_x((4, 4), (4, 2), rcond=1e-15)
# Test more unknowns than data
self.check_x((3, 10), (3, ), rcond=1e-15)
self.check_x((3, 10), (3, 3), rcond=1e-15)
# Test large rcond
self.check_x((4, 4), (4,), rcond=0.3)
self.check_x((2, 5), (2,), rcond=0.5)
self.check_x((5, 3), (5,), rcond=0.6)

def test_invalid_shape(self):
# test ValueError('Axis dimension mismatch')
self.check_shape((4, 3), (3, ))
self.check_shape((4, 3), (2, 2))
self.check_shape((4, 3), (3, 2, 3))
self.check_shape((4, 3), (10, 3, 3))
# Note that invalid shapes of a are already tested with pinv


@unittest.skipUnless(
cuda.cusolver_enabled, 'Only cusolver in CUDA 8.0 is supported')
@testing.gpu
ProTip! Use n and p to navigate between commits in a pull request.
You can’t perform that action at this time.