* update source files (#17279)

* fix - Invalid type for args: NDArray-or-Symbol-or-Scalar

* fix - unix gpu: use cudaMalloc

* fix - Not use CUDA_CALL

* fix - not Free(handle)

* fix - bug in c_lapack_api.h

python/mxnet/ndarray/numpy/linalg.py

@@ -21,7 +21,79 @@
 from . import _op as _mx_nd_np
 from . import _internal as _npi
 
-__all__ = ['norm', 'svd', 'cholesky', 'inv', 'det', 'slogdet', 'solve', 'tensorinv', 'tensorsolve']
+__all__ = ['norm', 'svd', 'cholesky', 'inv', 'det', 'slogdet', 'solve', 'tensorinv', 'tensorsolve', 'pinv']
+
+
+def pinv(a, rcond=1e-15, hermitian=False):
+    r"""
+    Compute the (Moore-Penrose) pseudo-inverse of a matrix.
+
+    Calculate the generalized inverse of a matrix using its
+    singular-value decomposition (SVD) and including all
+    *large* singular values.
+
+    Parameters
+    ----------
+    a : (..., M, N) ndarray
+        Matrix or stack of matrices to be pseudo-inverted.
+    rcond : (...) {float or ndarray of float}, optional
+        Cutoff for small singular values.
+        Singular values less than or equal to
+        ``rcond * largest_singular_value`` are set to zero.
+        Broadcasts against the stack of matrices.
+    hermitian : bool, optional
+        If True, `a` is assumed to be Hermitian (symmetric if real-valued),
+        enabling a more efficient method for finding singular values.
+        Defaults to False.
+
+    Returns
+    -------
+    B : (..., N, M) ndarray
+        The pseudo-inverse of `a`. If `a` is a `matrix` instance, then so
+        is `B`.
+
+    Raises
+    ------
+    MXNetError
+        If the SVD computation does not converge.
+
+    Notes
+    -----
+    The pseudo-inverse of a matrix A, denoted :math:`A^+`, is
+    defined as: "the matrix that 'solves' [the least-squares problem]
+    :math:`Ax = b`," i.e., if :math:`\\bar{x}` is said solution, then
+    :math:`A^+` is that matrix such that :math:`\\bar{x} = A^+b`.
+
+    It can be shown that if :math:`Q_1 \\Sigma Q_2^T = A` is the singular
+    value decomposition of A, then
+    :math:`A^+ = Q_2 \\Sigma^+ Q_1^T`, where :math:`Q_{1,2}` are
+    orthogonal matrices, :math:`\\Sigma` is a diagonal matrix consisting
+    of A's so-called singular values, (followed, typically, by
+    zeros), and then :math:`\\Sigma^+` is simply the diagonal matrix
+    consisting of the reciprocals of A's singular values
+    (again, followed by zeros). [1]_
+
+    References
+    ----------
+    .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
+           FL, Academic Press, Inc., 1980, pp. 139-142.
+
+    Examples
+    --------
+    The following example checks that ``a * a+ * a == a`` and
+    ``a+ * a * a+ == a+``:
+    >>> a = np.random.randn(2, 3)
+    >>> pinv_a = np.linalg.pinv(a)
+    >>> (a - np.dot(a, np.dot(pinv_a, a))).sum()
+    array(0.)
+    >>> (pinv_a - np.dot(pinv_a, np.dot(a, pinv_a))).sum()
+    array(0.)
+    """
+    if hermitian is True:
+        raise NotImplementedError("hermitian is not supported yet...")
+    if _mx_nd_np._np.isscalar(rcond):
+        return _npi.pinv_scalar_rcond(a, rcond, hermitian)
+    return _npi.pinv(a, rcond, hermitian)
 
 
 def norm(x, ord=None, axis=None, keepdims=False):

python/mxnet/numpy/linalg.py

@@ -20,7 +20,75 @@
 from __future__ import absolute_import
 from ..ndarray import numpy as _mx_nd_np
 
-__all__ = ['norm', 'svd', 'cholesky', 'inv', 'det', 'slogdet', 'solve', 'tensorinv', 'tensorsolve']
+__all__ = ['norm', 'svd', 'cholesky', 'inv', 'det', 'slogdet', 'solve', 'tensorinv', 'tensorsolve', 'pinv']
+
+
+def pinv(a, rcond=1e-15, hermitian=False):
+    r"""
+    Compute the (Moore-Penrose) pseudo-inverse of a matrix.
+
+    Calculate the generalized inverse of a matrix using its
+    singular-value decomposition (SVD) and including all
+    *large* singular values.
+
+    Parameters
+    ----------
+    a : (..., M, N) ndarray
+        Matrix or stack of matrices to be pseudo-inverted.
+    rcond : (...) {float or ndarray of float}, optional
+        Cutoff for small singular values.
+        Singular values less than or equal to
+        ``rcond * largest_singular_value`` are set to zero.
+        Broadcasts against the stack of matrices.
+    hermitian : bool, optional
+        If True, `a` is assumed to be Hermitian (symmetric if real-valued),
+        enabling a more efficient method for finding singular values.
+        Defaults to False.
+
+    Returns
+    -------
+    B : (..., N, M) ndarray
+        The pseudo-inverse of `a`. If `a` is a `matrix` instance, then so
+        is `B`.
+
+    Raises
+    ------
+    MXNetError
+        If the SVD computation does not converge.
+
+    Notes
+    -----
+    The pseudo-inverse of a matrix A, denoted :math:`A^+`, is
+    defined as: "the matrix that 'solves' [the least-squares problem]
+    :math:`Ax = b`," i.e., if :math:`\\bar{x}` is said solution, then
+    :math:`A^+` is that matrix such that :math:`\\bar{x} = A^+b`.
+
+    It can be shown that if :math:`Q_1 \\Sigma Q_2^T = A` is the singular
+    value decomposition of A, then
+    :math:`A^+ = Q_2 \\Sigma^+ Q_1^T`, where :math:`Q_{1,2}` are
+    orthogonal matrices, :math:`\\Sigma` is a diagonal matrix consisting
+    of A's so-called singular values, (followed, typically, by
+    zeros), and then :math:`\\Sigma^+` is simply the diagonal matrix
+    consisting of the reciprocals of A's singular values
+    (again, followed by zeros). [1]_
+
+    References
+    ----------
+    .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
+           FL, Academic Press, Inc., 1980, pp. 139-142.
+
+    Examples
+    --------
+    The following example checks that ``a * a+ * a == a`` and
+    ``a+ * a * a+ == a+``:
+    >>> a = np.random.randn(2, 3)
+    >>> pinv_a = np.linalg.pinv(a)
+    >>> (a - np.dot(a, np.dot(pinv_a, a))).sum()
+    array(0.)
+    >>> (pinv_a - np.dot(pinv_a, np.dot(a, pinv_a))).sum()
+    array(0.)
+    """
+    return _mx_nd_np.linalg.pinv(a, rcond, hermitian)
 
 
 def norm(x, ord=None, axis=None, keepdims=False):

python/mxnet/numpy_dispatch_protocol.py

@@ -148,6 +148,7 @@ def _run_with_array_ufunc_proto(*args, **kwargs):
     'linalg.solve',
     'linalg.tensorinv',
     'linalg.tensorsolve',
+    'linalg.pinv',
     'shape',
     'trace',
     'tril',

python/mxnet/symbol/numpy/linalg.py

@@ -22,7 +22,78 @@
 from . import _op as _mx_sym_np
 from . import _internal as _npi
 
-__all__ = ['norm', 'svd', 'cholesky', 'inv', 'det', 'slogdet', 'solve', 'tensorinv', 'tensorsolve']
+__all__ = ['norm', 'svd', 'cholesky', 'inv', 'det', 'slogdet', 'solve', 'tensorinv', 'tensorsolve', 'pinv']
+
+def pinv(a, rcond=1e-15, hermitian=False):
+    r"""
+    Compute the (Moore-Penrose) pseudo-inverse of a matrix.
+
+    Calculate the generalized inverse of a matrix using its
+    singular-value decomposition (SVD) and including all
+    *large* singular values.
+
+    Parameters
+    ----------
+    a : (..., M, N) ndarray
+        Matrix or stack of matrices to be pseudo-inverted.
+    rcond : (...) {float or ndarray of float}, optional
+        Cutoff for small singular values.
+        Singular values less than or equal to
+        ``rcond * largest_singular_value`` are set to zero.
+        Broadcasts against the stack of matrices.
+    hermitian : bool, optional
+        If True, `a` is assumed to be Hermitian (symmetric if real-valued),
+        enabling a more efficient method for finding singular values.
+        Defaults to False.
+
+    Returns
+    -------
+    B : (..., N, M) ndarray
+        The pseudo-inverse of `a`. If `a` is a `matrix` instance, then so
+        is `B`.
+
+    Raises
+    ------
+    MXNetError
+        If the SVD computation does not converge.
+
+    Notes
+    -----
+    The pseudo-inverse of a matrix A, denoted :math:`A^+`, is
+    defined as: "the matrix that 'solves' [the least-squares problem]
+    :math:`Ax = b`," i.e., if :math:`\\bar{x}` is said solution, then
+    :math:`A^+` is that matrix such that :math:`\\bar{x} = A^+b`.
+
+    It can be shown that if :math:`Q_1 \\Sigma Q_2^T = A` is the singular
+    value decomposition of A, then
+    :math:`A^+ = Q_2 \\Sigma^+ Q_1^T`, where :math:`Q_{1,2}` are
+    orthogonal matrices, :math:`\\Sigma` is a diagonal matrix consisting
+    of A's so-called singular values, (followed, typically, by
+    zeros), and then :math:`\\Sigma^+` is simply the diagonal matrix
+    consisting of the reciprocals of A's singular values
+    (again, followed by zeros). [1]_
+
+    References
+    ----------
+    .. [1] G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando,
+           FL, Academic Press, Inc., 1980, pp. 139-142.
+
+    Examples
+    --------
+    The following example checks that ``a * a+ * a == a`` and
+    ``a+ * a * a+ == a+``:
+    >>> a = np.random.randn(2, 3)
+    >>> pinv_a = np.linalg.pinv(a)
+    >>> (a - np.dot(a, np.dot(pinv_a, a))).sum()
+    array(0.)
+    >>> (pinv_a - np.dot(pinv_a, np.dot(a, pinv_a))).sum()
+    array(0.)
+    """
+    if hermitian is True:
+        raise NotImplementedError("hermitian is not supported yet...")
+    if _symbol._np.isscalar(rcond):
+        return _npi.pinv_scalar_rcond(a, rcond, hermitian)
+    return _npi.pinv(a, rcond, hermitian)
 
 
 def norm(x, ord=None, axis=None, keepdims=False):

src/operator/c_lapack_api.cc

@@ -78,6 +78,16 @@
     return 1; \
   }
 
+  #define MXNET_LAPACK_CWRAPPER9(func, dtype) \
+  int MXNET_LAPACK_##func(int matrix_layout, int m, int n, \
+                          dtype *a, int lda, dtype *s, \
+                          dtype *u, int ldu, \
+                          dtype *vt, int ldvt, \
+                          dtype *work, int lwork, int *iwork) { \
+    LOG(FATAL) << "MXNet build without lapack. Function " << #func << " is not available."; \
+    return 1; \
+  }
+
   #define MXNET_LAPACK_UNAVAILABLE(func) \
   int mxnet_lapack_##func(...) { \
     LOG(FATAL) << "MXNet build without lapack. Function " << #func << " is not available."; \
@@ -111,4 +121,7 @@
   MXNET_LAPACK_CWRAPPER7(sgesv, float)
   MXNET_LAPACK_CWRAPPER7(dgesv, double)
 
+  MXNET_LAPACK_CWRAPPER9(sgesdd, float)
+  MXNET_LAPACK_CWRAPPER9(dgesdd, double)
+
 #endif  // MSHADOW_USE_MKL == 0

src/operator/c_lapack_api.h

@@ -163,6 +163,23 @@ extern "C" {
 
   MXNET_LAPACK_FSIG_GESV(sgesv, float)
   MXNET_LAPACK_FSIG_GESV(dgesv, double)
+
+  #ifdef __ANDROID__
+    #define MXNET_LAPACK_FSIG_GESDD(func, dtype) \
+    int func##_(char *jobz, int *m, int *n, dtype *a, int *lda, dtype *s, \
+                dtype *u, int *ldu, \
+                dtype *vt, int *ldvt, \
+                dtype *work, int *lwork, int *iwork, int *info);
+  #else
+    #define MXNET_LAPACK_FSIG_GESDD(func, dtype) \
+    void func##_(char *jobz, int *m, int *n, dtype *a, int *lda, dtype *s, \
+                 dtype *u, int *ldu, \
+                 dtype *vt, int *ldvt, \
+                 dtype *work, int *lwork, int *iwork, int *info);
+  #endif
+
+  MXNET_LAPACK_FSIG_GESDD(sgesdd, float)
+  MXNET_LAPACK_FSIG_GESDD(dgesdd, double)
 }
 
 #endif  // MSHADOW_USE_MKL == 0
@@ -282,6 +299,24 @@ inline void flip(int m, int n, DType *b, int ldb, DType *a, int lda) {
   MXNET_LAPACK_CWRAP_GESVD(s, float)
   MXNET_LAPACK_CWRAP_GESVD(d, double)
 
+  // Computes the singular value decomposition of a general rectangular matrix
+  // using a divide and conquer method.
+  #define MXNET_LAPACK_CWRAP_GESDD(prefix, dtype) \
+  inline int MXNET_LAPACK_##prefix##gesdd(int matrix_layout, int m, int n, \
+                                          dtype *a, int lda, dtype *s, \
+                                          dtype *u, int ldu, \
+                                          dtype *vt, int ldvt, \
+                                          dtype *work, int lwork, int *iwork) { \
+    if (lwork != -1) { \
+      return LAPACKE_##prefix##gesdd(matrix_layout, 'O', m, n, a, lda, \
+                                     s, u, ldu, vt, ldvt); \
+    } \
+    *work = 0; \
+    return 0; \
+  }
+  MXNET_LAPACK_CWRAP_GESDD(s, float)
+  MXNET_LAPACK_CWRAP_GESDD(d, double)
+
   #define MXNET_LAPACK_CWRAP_GETRI(prefix, dtype) \
   inline int MXNET_LAPACK_##prefix##getri(int matrix_layout, int n, dtype *a, int lda, \
                                           int *ipiv, dtype *work, int lwork) { \
@@ -421,6 +456,25 @@ inline void flip(int m, int n, DType *b, int ldb, DType *a, int lda) {
   MXNET_LAPACK_CWRAP_GESVD(sgesvd, float)
   MXNET_LAPACK_CWRAP_GESVD(dgesvd, double)
 
+  #define MXNET_LAPACK_CWRAP_GESDD(func, dtype) \
+  inline int MXNET_LAPACK_##func(int matrix_layout, int m, int n, \
+                                 dtype *a, int lda, dtype *s, \
+                                 dtype *u, int ldu, \
+                                 dtype *vt, int ldvt, \
+                                 dtype *work, int lwork, int *iwork) { \
+    if (matrix_layout == MXNET_LAPACK_ROW_MAJOR) { \
+      CHECK(false) << "MXNET_LAPACK_" << #func << " implemented for row-major layout only"; \
+      return 1; \
+    } else { \
+      int info(0); \
+      char jobz('O'); \
+      func##_(&jobz, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, iwork, &info); \
+      return info; \
+    } \
+  }
+  MXNET_LAPACK_CWRAP_GESDD(sgesdd, float)
+  MXNET_LAPACK_CWRAP_GESDD(dgesdd, double)
+
   #define MXNET_LAPACK
 
   // Note: Both MXNET_LAPACK_*getrf, MXNET_LAPACK_*getri can only be called with col-major format
@@ -506,6 +560,13 @@ inline void flip(int m, int n, DType *b, int ldb, DType *a, int lda) {
   int MXNET_LAPACK_##func(int matrix_order, int n, int nrhs, dtype *a, \
                           int lda, int *ipiv, dtype *b, int ldb); \
 
+  #define MXNET_LAPACK_CWRAPPER9(func, dtype) \
+  int MXNET_LAPACK_##func(int matrix_layout, int m, int n, \
+                          dtype *a, int lda, dtype *s, \
+                          dtype *u, int ldu, \
+                          dtype *vt, int ldvt, \
+                          dtype *work, int lwork, int *iwork);
+
   #define MXNET_LAPACK_UNAVAILABLE(func) \
   int mxnet_lapack_##func(...);
   MXNET_LAPACK_CWRAPPER1(spotrf, float)
@@ -536,6 +597,9 @@ inline void flip(int m, int n, DType *b, int ldb, DType *a, int lda) {
   MXNET_LAPACK_CWRAPPER7(sgesv, float)
   MXNET_LAPACK_CWRAPPER7(dgesv, double)
 
+  MXNET_LAPACK_CWRAPPER9(sgesdd, float)
+  MXNET_LAPACK_CWRAPPER9(dgesdd, double)
+
   #undef MXNET_LAPACK_CWRAPPER1
   #undef MXNET_LAPACK_CWRAPPER2
   #undef MXNET_LAPACK_CWRAPPER3