• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1# Copyright 2022 Huawei Technologies Co., Ltd
2#
3# Licensed under the Apache License, Version 2.0 (the "License");
4# you may not use this file except in compliance with the License.
5# You may obtain a copy of the License at
6#
7# http://www.apache.org/licenses/LICENSE-2.0
8#
9# Unless required by applicable law or agreed to in writing, software
10# distributed under the License is distributed on an "AS IS" BASIS,
11# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12# See the License for the specific language governing permissions and
13# limitations under the License.
14# ============================================================================
15
16"""Operators for linalg function."""
17from __future__ import absolute_import
18
19import mindspore.ops as ops
20from mindspore.common import dtype as mstype
21from mindspore.common.tensor import Tensor
22from mindspore.ops import operations as P
23from mindspore.ops import functional as F
24from mindspore.ops.operations import _inner_ops as inner
25from mindspore.ops.function.math_func import _check_input_dtype, _check_attr_dtype
26from mindspore._c_expression import Tensor as Tensor_
27from mindspore.ops.auto_generate import geqrf
28
29from ..operations import linalg_ops
30from .._primitive_cache import _get_cache_prim
31
32__all__ = ['cond', 'eig', 'eigvals', 'geqrf', 'svd', 'pinv', 'qr']
33
34dtype_ = P.DType()
35geqrf_ = P.Geqrf()
36slice_ = P.Slice()
37
38def cond(A, p=None):
39    r"""
40    Returns the matrix norm or vector norm of a given tensor.
41
42    `p` is the calculation mode of norm. The following norm modes are supported.
43
44    ==========================  ================================ ==========================================
45    `p`                         norm for matrices                 norm for vectors
46    ==========================  ================================ ==========================================
47    ``None`` (default)          `2`-norm (see below)              `2`-norm (see below)
48    ``'fro'``                   Frobenius norm                    -- not supported --
49    ``'nuc'``                   nuclear norm                      -- not supported --
50    ``inf``                     :math:`max(sum(abs(x), dim=1))`   :math:`max(abs(x))`
51    ``-inf``                    :math:`min(sum(abs(x), dim=1))`   :math:`min(abs(x))`
52    ``0``                       -- not supported --               :math:`sum(x != 0)`
53    ``1``                       :math:`max(sum(abs(x), dim=0))`   as below
54    ``-1``                      :math:`min(sum(abs(x), dim=0))`   as below
55    ``2``                       largest singular value            as below
56    ``-2``                      smallest singular value           as below
57    other ``int`` or ``float``  -- not supported --               :math:`sum(abs(x)^{p})^{(1 / p)}`
58    ==========================  ================================ ==========================================
59
60    Note:
61        Currently, complex numbers are not supported.
62
63    Args:
64        A (Tensor): Tensor of shape :math:`(*, n)` or :math:`(*, m, n)`
65            where :math:`*` is zero or more batch dimensions.
66        p (Union[int, float, inf, -inf, 'fro', 'nuc'], optional): norm's mode. Refer to the table above for
67            behavior. Default: ``None``.
68
69    Returns:
70        Tensor, the result of norm calculation on the specified dimension, `dim`, has the same dtype as `A`.
71
72    Raises:
73        TypeError: If `A` is a vector and `p` is a str.
74        ValueError: If `A` is a matrices and `p` is not in valid mode.
75        ValueError: If `A` is a matrix and `p` is an integer that is not in [1, -1, 2, -2].
76
77    Supported Platforms:
78        ``GPU`` ``CPU``
79
80    Examples:
81        >>> import mindspore as ms
82        >>> x = ms.Tensor([[1.0, 0.0, -1.0], [0.0, 1.0, 0.0], [1.0, 0.0, 1.0]])
83        >>> print(ms.ops.cond(x))
84        1.4142
85        >>> print(ms.ops.cond(x, 'fro'))
86        3.1622777
87    """
88    matrix_inverse = _get_cache_prim(P.MatrixInverse)(adjoint=False)
89    if p is None:
90        p = 2
91    if A.dim() >= 3:
92        shape_ori = A.shape[0:-2]
93        A_flatten = ops.flatten(A, start_dim=0, end_dim=-3)
94        out = []
95        for i in range(A_flatten.shape[0]):
96            norm_a = F.norm(A_flatten[i], p)
97            norm_inv_a = F.norm(matrix_inverse(A_flatten[i]), p)
98            cond_i = ops.fill(mstype.float32, (1, 1), norm_a * norm_inv_a)
99            out.append(cond_i)
100        out_stacked = ops.hstack(out)
101        output = ops.reshape(out_stacked, shape_ori)
102        return output
103    norm_a = F.norm(A, p)
104    norm_inv_a = F.norm(matrix_inverse(A), p)
105    return norm_a * norm_inv_a
106
107
108def eig(A):
109    """
110    Computes the eigenvalues and eigenvectors of a square matrix(batch square matrices).
111
112    .. warning::
113        This is an experimental API that is subject to change or deletion.
114
115    Args:
116        A (Tensor): Square matrices of shape :math:`(*, N, N)`,
117            with float32, float64, complex64 or complex128 data type.
118
119    Returns:
120        - **eigen_values** (Tensor) - Shape :math:`(*, N)`. eigenvalues of
121          the corresponding matrix. The eigenvalues may not have an order.
122        - **eigen_vectors** (Tensor) - Shape :math:`(*, N, N)`,columns of eigen vectors represent
123        - **normalized** (unit length) eigenvectors of corresponding eigenvalues.
124
125    Raises:
126        TypeError: If dtype of `A` is not one of: float64, float32, complex64 or complex128.
127        TypeError: If `A` is not a Tensor.
128        ValueError: If `A` is not a square(batch squares).
129
130    Supported Platforms:
131        ``Ascend`` ``CPU``
132
133    Examples:
134        >>> import mindspore
135        >>> import numpy as np
136        >>> from mindspore import Tensor, ops
137        >>> input_x = Tensor(np.array([[1.0, 0.0], [0.0, 2.0]]), mindspore.float32)
138        >>> u, v = ops.eig(input_x)
139        >>> print(u)
140        [1.+0.j 2.+0.j]
141        >>> print(v)
142        [[1.+0.j 0.+0.j]
143         [0.+0.j 1.+0.j]]
144    """
145    return _get_cache_prim(P.Eig)(compute_v=True)(A)
146
147
148def eigvals(A):
149    """
150    Computes the eigenvalues of a square matrix(batch square matrices).
151
152    .. warning::
153        This is an experimental API that is subject to change or deletion.
154
155    Args:
156        A (Tensor): Square matrices of shape :math:`(*, N, N)`,
157            with float32, float64, complex64 or complex128 data type.
158
159    Returns:
160        Tensor, with shape :math:`(*, N)`. Returns the eigenvalues of
161        the corresponding matrix, which may not have an order.
162
163    Raises:
164        TypeError: If dtype of `A` is not one of: float64, float32, complex64 or complex128.
165        TypeError: If `A` is not a Tensor.
166        ValueError: If `A` is not a square(batch squares).
167
168    Supported Platforms:
169        ``Ascend`` ``CPU``
170
171    Examples:
172        >>> import mindspore
173        >>> from mindspore import Tensor, ops
174        >>> import numpy as np
175        >>> input_x = Tensor(np.array([[1.0, 0.0], [0.0, 2.0]]), mindspore.float32)
176        >>> u = ops.eigvals(input_x)
177        >>> print(u)
178        [1.+0.j 2.+0.j]
179    """
180    u, _ = _get_cache_prim(P.Eig)(compute_v=False)(A)
181    return u
182
183
184def svd(input, full_matrices=False, compute_uv=True):
185    """
186    Computes the singular value decompositions of one or more matrices.
187
188    If :math:`A` is a matrix, the svd returns the singular values :math:`S`, the left singular vectors :math:`U`
189    and the right singular vectors :math:`V`. It meets:
190
191    .. math::
192        A=U*diag(S)*V^{T}
193
194    Args:
195        input (Tensor): Tensor of the matrices to be decomposed. The shape should be :math:`(*, M, N)`,
196          the supported dtype are float32 and float64.
197        full_matrices (bool, optional): If true, compute full-sized :math:`U` and :math:`V`. If false, compute
198                                        only the leading P singular vectors, with P is the minimum of M and N.
199                                        Default: ``False`` .
200        compute_uv (bool, optional): If true, compute the left and right singular vectors.
201                                     If false, compute only the singular values. Default: ``True`` .
202
203    Returns:
204        - **s**  (Tensor) - Singular values. The shape is :math:`(*, P)`.
205        - **u**  (Tensor) - Left singular vectors. If `compute_uv` is False, u will not be returned.
206          The shape is :math:`(*, M, P)`. If `full_matrices` is True, the shape will be :math:`(*, M, M)`.
207        - **v**  (Tensor) - Right singular vectors. If `compute_uv` is False, v will not be returned.
208          The shape is :math:`(*, N, P)`. If `full_matrices` is True, the shape will be :math:`(*, N, N)`.
209
210    Raises:
211        TypeError: If `full_matrices` or `compute_uv` is not the type of bool.
212        TypeError: If the rank of input less than 2.
213        TypeError: If the type of input is not one of the following dtype: float32, float64.
214
215    Supported Platforms:
216        ``GPU`` ``CPU``
217
218    Examples:
219        >>> import numpy as np
220        >>> from mindspore import Tensor, set_context
221        >>> from mindspore import ops
222        >>> set_context(device_target="CPU")
223        >>> input = Tensor(np.array([[1, 2], [-4, -5], [2, 1]]).astype(np.float32))
224        >>> s, u, v = ops.svd(input, full_matrices=True, compute_uv=True)
225        >>> print(s)
226        [7.0652843 1.040081 ]
227        >>> print(u)
228        [[ 0.30821905 -0.48819482 0.81649697]
229         [-0.90613353  0.11070572 0.40824813]
230         [ 0.2896955   0.8656849  0.4082479 ]]
231        >>> print(v)
232        [[ 0.63863593 0.769509  ]
233         [ 0.769509  -0.63863593]]
234    """
235    svd_ = _get_cache_prim(linalg_ops.Svd)(full_matrices=full_matrices, compute_uv=compute_uv)
236
237    if compute_uv:
238        return svd_(input)
239
240    s, _, _ = svd_(input)
241    return s
242
243
244def pinv(x, *, atol=None, rtol=None, hermitian=False):
245    r"""
246    Computes the (Moore-Penrose) pseudo-inverse of a matrix.
247    This function is computed using SVD. If :math:`x=U*S*V^{T}` ,Than the pseudo-inverse of x is:
248    :math:`x^{+}=V*S^{+}*U^{T}` , :math:`S^{+}` is the reciprocal of each non-zero element on
249    the diagonal of S, and zero remains in place.
250    Batch matrices are supported. If x is a batch matrix, the output has the same batch dimension when
251    atol or rtol is float.
252    If atol or rtol is a Tensor, its shape must be broadcast to the singular value returned by
253    `x.svd <https://www.mindspore.cn/docs/en/master/api_python/ops/mindspore.ops.svd.html>`_ .
254    If x.shape is :math:`(B, M, N)`, and the shape of atol or rtol is :math:`(K, B)`, the output
255    shape is :math:`(K, B, N, M)`.
256    When the Hermitian is True, temporary support only real domain, x is treated as a real symmetric, so x is
257    not checked internally, and only use the lower triangular part in the computations.
258    When the singular value of x (or the norm of the eigenvalues when hermitian = True) that are below
259    threshold (:math:`max(atol, \sigma \cdot rtol)`, :math:`\sigma` as the largest singular value or
260    characteristic value), it is set to zero, and is not used in the computations.
261    If rtol is not specified and x is a matrix of dimensions (M, N), then rtol is set to
262    be :math:`rtol=max(M, N)*\varepsilon`, :math:`\varepsilon` is the
263    `eps <https://www.mindspore.cn/docs/en/master/api_python/ops/mindspore.ops.Eps.html>`_ value of x.dtype.
264    If rtol is not specified and atol specifies a value larger than zero, rtol is set to zero.
265
266    .. note::
267        This function uses
268        `svd <https://www.mindspore.cn/docs/en/master/api_python/ops/mindspore.ops.svd.html>`_ internally,
269        (or `eigh <https://www.mindspore.cn/docs/en/master/api_python/scipy/mindspore.scipy.linalg.eigh.html>`_ ,
270        when `hermitian = True` ). So it has the same problem as these functions. For details,
271        see the warnings in svd() and eigh().
272
273    Args:
274        x (Tensor): A matrix to be calculated. Only `float32`, `float64` are supported Tensor dtypes.
275            shape is :math:`(*, M, N)`, * is zero or more batch dimensions.
276
277            - When `hermitian` is ``True``, batch dimensions are not supported temporarily.
278
279    Keyword args:
280        atol (float, Tensor): absolute tolerance value. Default: ``None`` .
281        rtol (float, Tensor): relative tolerance value. Default: ``None`` .
282        hermitian (bool): An optional bool. x is assumed to be symmetric if real. Default: ``False`` .
283
284    Outputs:
285        - **output** (Tensor) - same type as input. Shape is :math:`(*, N, M)`, * is zero or more batch dimensions.
286
287    Raises:
288        TypeError: If `hermitian` is not a bool.
289        TypeError: If `x` is not a Tensor.
290        ValueError: If the dimension of `x` is less than 2.
291
292    Supported Platforms:
293        ``CPU``
294
295    Examples:
296        >>> import mindspore
297        >>> from mindspore import Tensor, ops
298        >>> x = Tensor([[4., 0.], [0., 5.]], mindspore.float32)
299        >>> output = ops.pinv(x)
300        >>> print(output)
301        [[0.25 0.  ]
302         [0.   0.2 ]]
303    """
304    if not isinstance(x, (Tensor, Tensor_)):
305        raise TypeError("The input x must be tensor")
306    if x.shape == ():
307        raise TypeError("For pinv, the 0-D input is not supported")
308    x_shape = F.shape(x)
309    if len(x_shape) < 2:
310        raise ValueError("input x should have 2 or more dimensions, " f"but got {len(x_shape)}.")
311    x_dtype = dtype_(x)
312    _check_input_dtype("x", x_dtype, [mstype.float32, mstype.float64], "pinv")
313    _check_attr_dtype("hermitian", hermitian, [bool], "pinv")
314
315    if atol is not None:
316        if rtol is None:
317            rtol = Tensor(0.0)
318    else:
319        atol = Tensor(0.0)
320        if rtol is None:
321            rtol = max(ops.Shape()(x)) * ops.Eps()(Tensor(1.0, x_dtype))
322
323    if not inner.IsInstance()(rtol, mstype.tensor_type):
324        rtol = F.cast(rtol, mstype.float32)
325    if not inner.IsInstance()(atol, mstype.tensor_type):
326        atol = F.cast(atol, mstype.float32)
327
328    if not hermitian:
329        s, u, v = linalg_ops.Svd()(x)
330        max_singular_val = _narrow(s, -1, 0, 1)
331        threshold = ops.Maximum()(atol.expand_dims(-1), rtol.expand_dims(-1) * max_singular_val)
332        condition = s > threshold
333        reciprocal_s_before = ops.Reciprocal()(s).broadcast_to(condition.shape)
334        zero = F.zeros(condition.shape, s.dtype)
335        s_pseudoinv = ops.Select()(condition, reciprocal_s_before, zero)
336        output = ops.matmul(v * s_pseudoinv.expand_dims(-2), _nd_transpose(ops.Conj()(u)))
337    else:
338        s, u = linalg_ops.Eigh()(x)
339        s_abs = ops.Abs()(s)
340        max_singular_val = ops.amax(s_abs, -1, True)
341        threshold = ops.Maximum()(atol.expand_dims(-1), rtol.expand_dims(-1) * max_singular_val)
342        condition = s_abs > threshold
343        reciprocal_s_before = ops.Reciprocal()(s)
344        zero = F.zeros(condition.shape, s.dtype)
345        s_pseudoinv = ops.Select()(condition, reciprocal_s_before, zero)
346        output = ops.matmul(u * s_pseudoinv.expand_dims(-2), _nd_transpose(ops.Conj()(u)))
347    return output
348
349
350def qr(input, mode='reduced'):
351    """
352    Returns the QR decomposition of one or more matrices.
353    If `mode` is 'reduced'(the default), compute the P columns of Q where P is minimum of the 2 innermost dimensions of
354    input. If `mode` is 'complete', compute full-sized Q and R.
355
356    Args:
357        input (Tensor): A matrix to be calculated. The matrix must be at least two dimensions, the supported dtype are
358            float16, float32, float64, complex64 and complex128.
359            Define the shape of input as :math:`(..., m, n)`, p as the
360            minimum values of m and n.
361        mode (Union['reduced', 'complete'], optional): If `mode` is ``'reduced'`` , computing reduce-sized QR
362            decomposition, otherwise, computing the full-sized QR decomposition. Default: ``'reduced'`` .
363
364    Returns:
365        - **Q** (Tensor) - The orthonormal matrices of input. If `mode` is 'complete', the shape is :math:`(m, m)`,
366          else the shape is :math:`(m, p)`. The dtype of `Q` is same as `input`.
367        - **R** (Tensor) - The upper triangular matrices of input. If `mode` is 'complete', the shape is :math:`(m, n)`,
368          else the shape is :math:`(p, n)`. The dtype of `R` is same as `input`.
369
370    Raises:
371        TypeError: If `input` is not a Tensor.
372        TypeError: If `mode` is neither 'reduced' nor 'complete'.
373        ValueError: If the dimension of `input` is less than 2.
374
375    Supported Platforms:
376        ``Ascend`` ``GPU`` ``CPU``
377
378    Examples:
379        >>> input = Tensor(np.array([[20., -31, 7], [4, 270, -90], [-8, 17, -32]]), mstype.float32)
380        >>> Q, R = ops.qr(input)
381        >>> print(Q)
382        [[-0.912871    0.16366126  0.37400758]
383         [-0.18257418 -0.9830709  -0.01544376]
384         [ 0.36514837 -0.08238228  0.92729706]]
385        >>> print(R)
386        [[ -21.908903  -14.788506  -1.6431675]
387        [    0.       -271.9031    92.25824  ]
388        [    0.          0.       -25.665514 ]]
389    """
390    if mode not in ('reduced', 'complete'):
391        raise TypeError(f"For qr, the arg mode must be 'reduced' or 'complete', but got {mode}.")
392    qr_ = _get_cache_prim(P.Qr)(mode == 'complete')
393    return qr_(input)
394
395
396def _narrow(x, axis, start, length):
397    begins = [0] * x.ndim
398    begins[axis] = start
399    sizes = list(x.shape)
400    sizes[axis] = length
401    return slice_(x, begins, sizes)
402
403
404def _nd_transpose(a):
405    """
406    _nd_transpose
407    """
408    dims = a.ndim
409    if dims < 2:
410        raise TypeError("to do _nd_transpose for input a's ndim is not greater or equal to 2d, which is invalid.")
411    axes = ops.make_range(0, dims)
412    axes = axes[:-2] + (axes[-1],) + (axes[-2],)
413    return ops.transpose(a, axes)
414
415
416__all__.sort()
417