diff --git a/python/mxnet/ndarray/numpy/linalg.py b/python/mxnet/ndarray/numpy/linalg.py index d443d68cda0d..74ba41f22979 100644 --- a/python/mxnet/ndarray/numpy/linalg.py +++ b/python/mxnet/ndarray/numpy/linalg.py @@ -21,7 +21,7 @@ from . import _op as _mx_nd_np from . import _internal as _npi -__all__ = ['norm', 'svd', 'cholesky', 'inv'] +__all__ = ['norm', 'svd', 'cholesky', 'inv', 'det', 'slogdet'] def norm(x, ord=None, axis=None, keepdims=False): @@ -242,3 +242,113 @@ def inv(a): [ 0.75000006, -0.25000003]]]) """ return _npi.inv(a) + + +def det(a): + r""" + Compute the determinant of an array. + + Parameters + ---------- + a : (..., M, M) ndarray + Input array to compute determinants for. + + Returns + ------- + det : (...) ndarray + Determinant of `a`. + + See Also + -------- + slogdet : Another way to represent the determinant, more suitable + for large matrices where underflow/overflow may occur. + + Notes + ----- + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. + The determinant is computed via LU factorization using the LAPACK + routine z/dgetrf. + + Examples + -------- + The determinant of a 2-D array [[a, b], [c, d]] is ad - bc: + >>> a = np.array([[1, 2], [3, 4]]) + >>> np.linalg.det(a) + -2.0 + + Computing determinants for a stack of matrices: + >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ]) + >>> a.shape + (3, 2, 2) + + >>> np.linalg.det(a) + array([-2., -3., -8.]) + """ + return _npi.det(a) + + +def slogdet(a): + r""" + Compute the sign and (natural) logarithm of the determinant of an array. + If an array has a very small or very large determinant, then a call to + `det` may overflow or underflow. This routine is more robust against such + issues, because it computes the logarithm of the determinant rather than + the determinant itself. + + Parameters + ---------- + a : (..., M, M) ndarray + Input array, has to be a square 2-D array. + + Returns + ------- + sign : (...) ndarray + A number representing the sign of the determinant. For a real matrix, + this is 1, 0, or -1. + logdet : (...) array_like + The natural log of the absolute value of the determinant. + If the determinant is zero, then `sign` will be 0 and `logdet` will be + -Inf. In all cases, the determinant is equal to ``sign * np.exp(logdet)``. + + See Also + -------- + det + + Notes + ----- + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. + The determinant is computed via LU factorization using the LAPACK + routine z/dgetrf. + + Examples + -------- + The determinant of a 2-D array ``[[a, b], [c, d]]`` is ``ad - bc``: + >>> a = np.array([[1, 2], [3, 4]]) + >>> (sign, logdet) = np.linalg.slogdet(a) + >>> (sign, logdet) + (-1., 0.69314718055994529) + + >>> sign * np.exp(logdet) + -2.0 + + Computing log-determinants for a stack of matrices: + >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ]) + >>> a.shape + (3, 2, 2) + + >>> sign, logdet = np.linalg.slogdet(a) + >>> (sign, logdet) + (array([-1., -1., -1.]), array([ 0.69314718, 1.09861229, 2.07944154])) + + >>> sign * np.exp(logdet) + array([-2., -3., -8.]) + + This routine succeeds where ordinary `det` does not: + >>> np.linalg.det(np.eye(500) * 0.1) + 0.0 + >>> np.linalg.slogdet(np.eye(500) * 0.1) + (1., -1151.2925464970228) + """ + return _npi.slogdet(a) diff --git a/python/mxnet/numpy/linalg.py b/python/mxnet/numpy/linalg.py index b0552ee9f319..fbe3631eb6e6 100644 --- a/python/mxnet/numpy/linalg.py +++ b/python/mxnet/numpy/linalg.py @@ -20,7 +20,7 @@ from __future__ import absolute_import from ..ndarray import numpy as _mx_nd_np -__all__ = ['norm', 'svd', 'cholesky', 'inv'] +__all__ = ['norm', 'svd', 'cholesky', 'inv', 'det', 'slogdet'] def norm(x, ord=None, axis=None, keepdims=False): @@ -260,3 +260,113 @@ def inv(a): [ 0.75000006, -0.25000003]]]) """ return _mx_nd_np.linalg.inv(a) + + +def det(a): + r""" + Compute the determinant of an array. + + Parameters + ---------- + a : (..., M, M) ndarray + Input array to compute determinants for. + + Returns + ------- + det : (...) ndarray + Determinant of `a`. + + See Also + -------- + slogdet : Another way to represent the determinant, more suitable + for large matrices where underflow/overflow may occur. + + Notes + ----- + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. + The determinant is computed via LU factorization using the LAPACK + routine z/dgetrf. + + Examples + -------- + The determinant of a 2-D array [[a, b], [c, d]] is ad - bc: + >>> a = np.array([[1, 2], [3, 4]]) + >>> np.linalg.det(a) + -2.0 + + Computing determinants for a stack of matrices: + >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ]) + >>> a.shape + (3, 2, 2) + + >>> np.linalg.det(a) + array([-2., -3., -8.]) + """ + return _mx_nd_np.linalg.det(a) + + +def slogdet(a): + r""" + Compute the sign and (natural) logarithm of the determinant of an array. + If an array has a very small or very large determinant, then a call to + `det` may overflow or underflow. This routine is more robust against such + issues, because it computes the logarithm of the determinant rather than + the determinant itself. + + Parameters + ---------- + a : (..., M, M) ndarray + Input array, has to be a square 2-D array. + + Returns + ------- + sign : (...) ndarray + A number representing the sign of the determinant. For a real matrix, + this is 1, 0, or -1. + logdet : (...) array_like + The natural log of the absolute value of the determinant. + If the determinant is zero, then `sign` will be 0 and `logdet` will be + -Inf. In all cases, the determinant is equal to ``sign * np.exp(logdet)``. + + See Also + -------- + det + + Notes + ----- + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. + The determinant is computed via LU factorization using the LAPACK + routine z/dgetrf. + + Examples + -------- + The determinant of a 2-D array ``[[a, b], [c, d]]`` is ``ad - bc``: + >>> a = np.array([[1, 2], [3, 4]]) + >>> (sign, logdet) = np.linalg.slogdet(a) + >>> (sign, logdet) + (-1., 0.69314718055994529) + + >>> sign * np.exp(logdet) + -2.0 + + Computing log-determinants for a stack of matrices: + >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ]) + >>> a.shape + (3, 2, 2) + + >>> sign, logdet = np.linalg.slogdet(a) + >>> (sign, logdet) + (array([-1., -1., -1.]), array([ 0.69314718, 1.09861229, 2.07944154])) + + >>> sign * np.exp(logdet) + array([-2., -3., -8.]) + + This routine succeeds where ordinary `det` does not: + >>> np.linalg.det(np.eye(500) * 0.1) + 0.0 + >>> np.linalg.slogdet(np.eye(500) * 0.1) + (1., -1151.2925464970228) + """ + return _mx_nd_np.linalg.slogdet(a) diff --git a/python/mxnet/symbol/numpy/linalg.py b/python/mxnet/symbol/numpy/linalg.py index 3c8118550559..cf33777b2637 100644 --- a/python/mxnet/symbol/numpy/linalg.py +++ b/python/mxnet/symbol/numpy/linalg.py @@ -22,7 +22,7 @@ from . import _op as _mx_sym_np from . import _internal as _npi -__all__ = ['norm', 'svd', 'cholesky', 'inv'] +__all__ = ['norm', 'svd', 'cholesky', 'inv', 'det', 'slogdet'] def norm(x, ord=None, axis=None, keepdims=False): @@ -229,3 +229,113 @@ def inv(a): [ 0.75000006, -0.25000003]]]) """ return _npi.inv(a) + + +def det(a): + r""" + Compute the determinant of an array. + + Parameters + ---------- + a : (..., M, M) ndarray + Input array to compute determinants for. + + Returns + ------- + det : (...) ndarray + Determinant of `a`. + + See Also + -------- + slogdet : Another way to represent the determinant, more suitable + for large matrices where underflow/overflow may occur. + + Notes + ----- + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. + The determinant is computed via LU factorization using the LAPACK + routine z/dgetrf. + + Examples + -------- + The determinant of a 2-D array [[a, b], [c, d]] is ad - bc: + >>> a = np.array([[1, 2], [3, 4]]) + >>> np.linalg.det(a) + -2.0 + + Computing determinants for a stack of matrices: + >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ]) + >>> a.shape + (3, 2, 2) + + >>> np.linalg.det(a) + array([-2., -3., -8.]) + """ + return _npi.det(a) + + +def slogdet(a): + r""" + Compute the sign and (natural) logarithm of the determinant of an array. + If an array has a very small or very large determinant, then a call to + `det` may overflow or underflow. This routine is more robust against such + issues, because it computes the logarithm of the determinant rather than + the determinant itself. + + Parameters + ---------- + a : (..., M, M) ndarray + Input array, has to be a square 2-D array. + + Returns + ------- + sign : (...) ndarray + A number representing the sign of the determinant. For a real matrix, + this is 1, 0, or -1. + logdet : (...) array_like + The natural log of the absolute value of the determinant. + If the determinant is zero, then `sign` will be 0 and `logdet` will be + -Inf. In all cases, the determinant is equal to ``sign * np.exp(logdet)``. + + See Also + -------- + det + + Notes + ----- + Broadcasting rules apply, see the `numpy.linalg` documentation for + details. + The determinant is computed via LU factorization using the LAPACK + routine z/dgetrf. + + Examples + -------- + The determinant of a 2-D array ``[[a, b], [c, d]]`` is ``ad - bc``: + >>> a = np.array([[1, 2], [3, 4]]) + >>> (sign, logdet) = np.linalg.slogdet(a) + >>> (sign, logdet) + (-1., 0.69314718055994529) + + >>> sign * np.exp(logdet) + -2.0 + + Computing log-determinants for a stack of matrices: + >>> a = np.array([ [[1, 2], [3, 4]], [[1, 2], [2, 1]], [[1, 3], [3, 1]] ]) + >>> a.shape + (3, 2, 2) + + >>> sign, logdet = np.linalg.slogdet(a) + >>> (sign, logdet) + (array([-1., -1., -1.]), array([ 0.69314718, 1.09861229, 2.07944154])) + + >>> sign * np.exp(logdet) + array([-2., -3., -8.]) + + This routine succeeds where ordinary `det` does not: + >>> np.linalg.det(np.eye(500) * 0.1) + 0.0 + >>> np.linalg.slogdet(np.eye(500) * 0.1) + (1., -1151.2925464970228) + """ + return _npi.slogdet(a) diff --git a/src/operator/tensor/la_op-inl.h b/src/operator/tensor/la_op-inl.h index 35edbbd42565..d580cced4ec5 100644 --- a/src/operator/tensor/la_op-inl.h +++ b/src/operator/tensor/la_op-inl.h @@ -502,6 +502,9 @@ struct det { static void op(const Tensor& A, const Tensor& det, const Tensor& LU, const Tensor& pivot, const OpContext& ctx, const nnvm::NodeAttrs& attrs) { + if (A.shape_.Size() == 0U) { + return; + } Stream *s = ctx.get_stream(); Tensor sign = ctx.requested[0] .get_space_typed(det.shape_, s); @@ -524,6 +527,9 @@ struct slogdet { const Tensor& logabsdet, const Tensor& LU, const Tensor& pivot, const OpContext& ctx, const nnvm::NodeAttrs& attrs) { + if (A.shape_.Size() == 0U) { + return; + } Stream *s = ctx.get_stream(); Copy(LU, A, s); linalg_batch_getrf(LU, pivot, false, s); @@ -921,6 +927,9 @@ struct det_backward { using namespace mshadow; using namespace mshadow::expr; using namespace mxnet_op; + if (dA.shape_.Size() == 0U) { + return; + } // compute inverse(A) and stores it to LU linalg_batch_det_backward_helper(LU, pivot, det, dA, DType(0), ctx); const_cast&>(dA) = broadcast_to(reshape(det * ddet, \ @@ -949,6 +958,9 @@ struct slogdet_backward { using namespace mshadow; using namespace mshadow::expr; using namespace mxnet_op; + if (dA.shape_.Size() == 0U) { + return; + } // compute inverse(A) and stores it to LU linalg_batch_det_backward_helper(LU, pivot, logabsdet, dA, DType(-INFINITY), ctx); const_cast&>(dA) = broadcast_to(reshape(dlogabsdet, \ diff --git a/src/operator/tensor/la_op.cc b/src/operator/tensor/la_op.cc index 8f9da23711a0..8407307bdd6d 100644 --- a/src/operator/tensor/la_op.cc +++ b/src/operator/tensor/la_op.cc @@ -945,6 +945,7 @@ NNVM_REGISTER_OP(_backward_linalg_inverse) NNVM_REGISTER_OP(_linalg_det) .add_alias("linalg_det") +.add_alias("_npi_det") .describe(R"code(Compute the determinant of a matrix. Input is a tensor *A* of dimension *n >= 2*. @@ -997,6 +998,7 @@ NNVM_REGISTER_OP(_backward_linalg_det) NNVM_REGISTER_OP(_linalg_slogdet) .add_alias("linalg_slogdet") +.add_alias("_npi_slogdet") .describe(R"code(Compute the sign and log of the determinant of a matrix. Input is a tensor *A* of dimension *n >= 2*. diff --git a/src/operator/tensor/la_op.h b/src/operator/tensor/la_op.h index e024693e3819..5fe7a92e2a12 100644 --- a/src/operator/tensor/la_op.h +++ b/src/operator/tensor/la_op.h @@ -26,6 +26,7 @@ #define MXNET_OPERATOR_TENSOR_LA_OP_H_ #include +#include #include #include #include "../mshadow_op.h" @@ -428,7 +429,11 @@ inline bool DetShape(const nnvm::NodeAttrs& attrs, CHECK_EQ(in[ndim-2], in[ndim-1]) << "Input A's last two dimension must be equal"; mxnet::TShape out; if (ndim == 2) { - out = mxnet::TShape(1, 1); + if (Imperative::Get()->is_np_shape() || in.Size() == 0U) { + out = mxnet::TShape(0, 1); + } else { + out = mxnet::TShape(1, 1); + } } else { out = mxnet::TShape(in.begin(), in.end() - 2); } @@ -896,19 +901,22 @@ void LaOpDetBackward(const nnvm::NodeAttrs& attrs, const std::vector& req, const std::vector& outputs) { using namespace mshadow; + if (outputs[0].shape_.Size() == 0U) { + return; + } Stream *s = ctx.get_stream(); CHECK_EQ(inputs.size(), onum + 3); CHECK_EQ(outputs.size(), 1); MSHADOW_SGL_DBL_TYPE_SWITCH(outputs[0].type_flag_, OType, { std::vector tspace(outputs); - for ( int i = 0; i < onum; ++i ) { + for ( size_t i = 0; i < outputs.size(); ++i ) { if ( req[i] == kAddTo ) { tspace[i].dptr_ = ctx.requested[0] .get_space_typed(Shape1(outputs[i].Size()), s).dptr_; } } LaOpDetBackwardCaller::op(inputs, tspace, attrs, ctx); - for ( int i = 0; i < onum; ++i ) { + for ( size_t i = 0; i < outputs.size(); ++i ) { if ( req[i] == kAddTo ) { Tensor out = outputs[i].FlatTo1D(s); out += tspace[i].FlatTo1D(s); diff --git a/tests/python/unittest/test_numpy_interoperability.py b/tests/python/unittest/test_numpy_interoperability.py index 87c48ab5fd79..206686f68de8 100644 --- a/tests/python/unittest/test_numpy_interoperability.py +++ b/tests/python/unittest/test_numpy_interoperability.py @@ -308,6 +308,16 @@ def _add_workload_linalg_inv(): OpArgMngr.add_workload('linalg.inv', np.array(_np.ones((0, 1, 1)), dtype=np.float64)) +def _add_workload_linalg_det(): + OpArgMngr.add_workload('linalg.det', np.array(_np.ones((2, 2)), dtype=np.float32)) + OpArgMngr.add_workload('linalg.det', np.array(_np.ones((0, 1, 1)), dtype=np.float64)) + + +def _add_workload_linalg_slogdet(): + OpArgMngr.add_workload('linalg.slogdet', np.array(_np.ones((2, 2)), dtype=np.float32)) + OpArgMngr.add_workload('linalg.slogdet', np.array(_np.ones((0, 1, 1)), dtype=np.float64)) + + def _add_workload_trace(): OpArgMngr.add_workload('trace', np.random.uniform(size=(4, 1))) OpArgMngr.add_workload('trace', np.random.uniform(size=(3, 2))) @@ -1332,6 +1342,8 @@ def _prepare_workloads(): _add_workload_linalg_norm() _add_workload_linalg_cholesky() _add_workload_linalg_inv() + _add_workload_linalg_det() + _add_workload_linalg_slogdet() _add_workload_trace() _add_workload_tril() _add_workload_outer() diff --git a/tests/python/unittest/test_numpy_op.py b/tests/python/unittest/test_numpy_op.py index dc99fc6fc251..f98783b13590 100644 --- a/tests/python/unittest/test_numpy_op.py +++ b/tests/python/unittest/test_numpy_op.py @@ -3373,6 +3373,106 @@ def check_inv(A_inv, data_np): check_inv(A_inv, data_np) +@with_seed() +@use_np +def test_np_linalg_det(): + class TestDet(HybridBlock): + def __init__(self): + super(TestDet, self).__init__() + + def hybrid_forward(self, F, a): + return F.np.linalg.det(a) + + # test non zero size input + tensor_shapes = [ + (2, 0, 2, 2), + (4, 4), + (0, 2, 2, 2), + (3, 3, 3), + (0, 2, 2), + (2, 2, 2, 2, 2), + (1, 1), + ] + types = [_np.float32, _np.float64] + grad_reqs = ['write', 'add', 'null'] + + for hybridize, dtype, shape, grad_req in itertools.product([True, False], types, tensor_shapes, grad_reqs): + a_shape = (1,) + shape + test_det = TestDet() + if hybridize: + test_det.hybridize() + a = rand_ndarray(shape=a_shape, dtype=dtype).as_np_ndarray() + a.attach_grad(grad_req) + np_out = _np.linalg.det(a.asnumpy()) + with mx.autograd.record(): + mx_out = test_det(a) + assert mx_out.shape == np_out.shape + assert_almost_equal(mx_out.asnumpy(), np_out, rtol=1e-1, atol=1e-1) + if grad_req != 'null': + print(shape, grad_req) + mx_out.backward() + + # Test imperative once again + mx_out = np.linalg.det(a) + np_out = _np.linalg.det(a.asnumpy()) + assert_almost_equal(mx_out.asnumpy(), np_out, rtol=1e-1, atol=1e-1) + + # test numeric gradient + a_sym = mx.sym.Variable("a").as_np_ndarray() + mx_sym = mx.sym.np.linalg.det(a_sym).as_nd_ndarray() + if 0 not in shape and grad_req != 'null': + check_numeric_gradient(mx_sym, [a.as_nd_ndarray()], + rtol=1e-1, atol=1e-1, dtype=dtype) + + +@with_seed() +@use_np +def test_np_linalg_slogdet(): + class TestSlogdet(HybridBlock): + def __init__(self): + super(TestSlogdet, self).__init__() + + def hybrid_forward(self, F, a): + return F.np.linalg.slogdet(a) + + # test non zero size input + tensor_shapes = [ + (2, 0, 2, 2), + (5, 5), + (0, 2, 2, 2), + (3, 3, 3), + (0, 3, 3), + (2, 2, 2, 2, 2), + (1, 1), + ] + + types = [_np.float32, _np.float64] + grad_reqs = ['write', 'add', 'null'] + + for hybridize, a_shape, dtype, grad_req in itertools.product([True, False], tensor_shapes, types, grad_reqs): + test_slogdet = TestSlogdet() + if hybridize: + test_slogdet.hybridize() + a = rand_ndarray(shape=a_shape, dtype=dtype).as_np_ndarray() + a.attach_grad(grad_req) + + np_out = _np.linalg.slogdet(a.asnumpy()) + with mx.autograd.record(): + mx_out = test_slogdet(a) + assert mx_out[0].shape == np_out[0].shape + assert mx_out[1].shape == np_out[1].shape + assert_almost_equal(mx_out[0].asnumpy(), np_out[0], rtol=1e-1, atol=1e-1) + assert_almost_equal(mx_out[1].asnumpy(), np_out[1], rtol=1e-1, atol=1e-1) + if grad_req != 'null': + mx_out[1].backward() + + # Test imperative once again + mx_out = np.linalg.slogdet(a) + np_out = _np.linalg.slogdet(a.asnumpy()) + assert_almost_equal(mx_out[0].asnumpy(), np_out[0], rtol=1e-1, atol=1e-1) + assert_almost_equal(mx_out[1].asnumpy(), np_out[1], rtol=1e-1, atol=1e-1) + + @with_seed() @use_np def test_np_vstack():