nlcpy.linalg.decomposition のソースコード

#
# * The source code in this file is based on the soure code of NumPy.
#
# # NLCPy License #
#
#     Copyright (c) 2020 NEC Corporation
#     All rights reserved.
#
#     Redistribution and use in source and binary forms, with or without
#     modification, are permitted provided that the following conditions are met:
#     * Redistributions of source code must retain the above copyright notice,
#       this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above copyright notice,
#       this list of conditions and the following disclaimer in the documentation
#       and/or other materials provided with the distribution.
#     * Neither NEC Corporation nor the names of its contributors may be
#       used to endorse or promote products derived from this software
#       without specific prior written permission.
#
#     THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
#     ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
#     WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
#     DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
#     FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
#     (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
#     LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
#     ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
#     (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
#     SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
# # NumPy License #
#
#     Copyright (c) 2005-2020, NumPy Developers.
#     All rights reserved.
#
#     Redistribution and use in source and binary forms, with or without
#     modification, are permitted provided that the following conditions are met:
#     * Redistributions of source code must retain the above copyright notice,
#       this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above copyright notice,
#       this list of conditions and the following disclaimer in the documentation
#       and/or other materials provided with the distribution.
#     * Neither the name of the NumPy Developers nor the names of any contributors may be
#       used to endorse or promote products derived from this software
#       without specific prior written permission.
#
#     THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
#     ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
#     WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
#     DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
#     FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
#     (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
#     LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
#     ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
#     (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
#     SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#

import nlcpy
import numpy
import warnings
from nlcpy import veo
from nlcpy.request import request
from . import util


def _take_along_axis(arr, indices, axis):
    shape = arr.shape
    shape_ones = (1,) * indices.ndim
    dest_dims = list(range(axis)) + [None] + list(range(axis + 1, indices.ndim))
    fancy_index = []
    for dim, n in zip(dest_dims, shape):
        if dim is None:
            fancy_index.append(indices)
        else:
            ind_shape = shape_ones[:dim] + (-1,) + shape_ones[dim + 1:]
            fancy_index.append(nlcpy.arange(n).reshape(ind_shape))
    return arr[tuple(fancy_index)]


[ドキュメント]def svd(a, full_matrices=True, compute_uv=True, hermitian=False): """Singular Value Decomposition. When *a* is a 2D array, it is factorized as ``u @ nlcpy.diag(s) @ vh = (u * s) @ vh``, where *u* and *vh* are 2D unitary arrays and *s* is a 1D array of *a*'s singular values. When *a* is higher-dimensional, SVD is applied in stacked mode as explained below. Parameters ---------- a : (..., M, N) array_like A real or complex array with a.ndim >= 2. full_matrices : bool, optional If True (default), *u* and *vh* have the shapes ``(..., M, M)`` and ``(..., N, N)``, respectively. Otherwise, the shapes are ``(..., M, K)`` and ``(..., K, N)``, respectively, where ``K = min(M, N)``. compute_uv : bool, optional Whether or not to compute *u* and *vh* in addition to *s*. True by default. 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 ------- u : {(..., M, M), (..., M, K)} ndarray Unitary array(s). The first ``a.ndim - 2`` dimensions have the same size as those of the input *a*. The size of the last two dimensions depends on the value of *full_matrices*. Only returned when *compute_uv* is True. s : (..., K) ndarray Vector(s) with the singular values, within each vector sorted in descending order. The first ``a.ndim - 2`` dimensions have the same size as those of the input *a*. vh : {(..., N, N), (..., K, N)} ndarray Unitary array(s). The first ``a.ndim - 2`` dimensions have the same size as those of the input *a*. The size of the last two dimensions depends on the value of *full_matrices*. Only returned when *compute_uv* is True. Note ---- The decomposition is performed using LAPACK routine ``_gesdd``. SVD is usually described for the factorization of a 2D matrix :math:`A`. The higher-dimensional case will be discussed below. In the 2D case, SVD is written as :math:`A=USV^{H}`, where :math:`A = a`, :math:`U = u`, :math:`S = nlcpy.diag(s)` and :math:`V^{H} = vh`. The 1D array `s` contains the singular values of `a` and `u` and `vh` are unitary. The rows of `vh` are the eigenvectors of :math:`A^{H}A` and the columns of `u` are the eigenvectors of :math:`AA^{H}`. In both cases the corresponding (possibly non-zero) eigenvalues are given by ``s**2``. If `a` has more than two dimensions, then broadcasting rules apply, as explained in :ref:`Linear algebra on several matrices at once <linalg_several_matrices_at_once>`. This means that SVD is working in "stacked" mode: it iterates over all indices of the first ``a.ndim - 2`` dimensions and for each combination SVD is applied to the last two indices. Examples -------- >>> import nlcpy as vp >>> from nlcpy import testing >>> a = vp.random.randn(9, 6) + 1j*vp.random.randn(9, 6) Reconstruction based on full SVD, 2D case: >>> u, s, vh = vp.linalg.svd(a, full_matrices=True) >>> u.shape, s.shape, vh.shape ((9, 9), (6,), (6, 6)) >>> vp.testing.assert_allclose(a, vp.dot(u[:, :6] * s, vh)) >>> smat = vp.zeros((9, 6), dtype=complex) >>> smat[:6, :6] = vp.diag(s) >>> vp.testing.assert_allclose(a, vp.dot(u, vp.dot(smat, vh))) Reconstruction based on reduced SVD, 2D case: >>> u, s, vh = vp.linalg.svd(a, full_matrices=False) >>> u.shape, s.shape, vh.shape ((9, 6), (6,), (6, 6)) >>> vp.testing.assert_allclose(a, vp.dot(u * s, vh)) >>> smat = vp.diag(s) >>> vp.testing.assert_allclose(a, vp.dot(u, vp.dot(smat, vh))) """ a = nlcpy.asarray(a) util._assertRankAtLeast2(a) if hermitian: util._assertNdSquareness(a) # used to match the contiguous of result to numpy. c_order = a.flags.c_contiguous or sum([i > 1 for i in a.shape[:-2]]) < 2 a_complex = a.dtype.char in 'FD' if a.dtype == 'F': dtype = 'F' f_dtype = 'f' elif a.dtype == 'D': dtype = 'D' f_dtype = 'd' elif a.dtype == 'f': dtype = 'f' f_dtype = 'f' else: dtype = 'd' f_dtype = 'd' if hermitian: if compute_uv: # lapack returns eigenvalues in reverse order, so to reconsist. s, u = nlcpy.linalg.eigh(a) signs = nlcpy.sign(s) s = abs(s) sidx = nlcpy.argsort(s)[..., ::-1] signs = _take_along_axis(signs, sidx, signs.ndim - 1) s = _take_along_axis(s, sidx, s.ndim - 1) u = _take_along_axis(u, sidx[..., None, :], u.ndim - 1) # singular values are unsigned, move the sign into v vt = nlcpy.conjugate(u * signs[..., None, :]) vt = nlcpy.moveaxis(vt, -2, -1) return u, s, vt else: s = nlcpy.linalg.eigvalsh(a) s = nlcpy.sort(abs(s))[..., ::-1] return s m = a.shape[-2] n = a.shape[-1] min_mn = min(m, n) max_mn = max(m, n) if a.size == 0: s = nlcpy.empty(a.shape[:-2] + (min_mn,), f_dtype) if compute_uv: if full_matrices: u_shape = a.shape[:-1] + (m, ) vt_shape = a.shape[:-2] + (n, n) else: u_shape = a.shape[:-1] + (min_mn, ) vt_shape = a.shape[:-2] + (min_mn, n) u = nlcpy.empty(u_shape, dtype=dtype) vt = nlcpy.empty(vt_shape, dtype=dtype) return u, s, vt else: return s a = nlcpy.array(nlcpy.moveaxis(a, (-1, -2), (1, 0)), dtype=dtype, order='F') if compute_uv: if full_matrices: u = nlcpy.empty((m, m) + a.shape[2:], dtype=dtype, order='F') vt = nlcpy.empty((n, n) + a.shape[2:], dtype=dtype, order='F') job = 'A' else: u = nlcpy.empty((m, m) + a.shape[2:], dtype=dtype, order='F') vt = nlcpy.empty((min_mn, n) + a.shape[2:], dtype=dtype, order='F') job = 'S' else: u = nlcpy.empty(1, dtype=dtype) vt = nlcpy.empty(1, dtype=dtype) job = 'N' if a_complex: mnthr1 = int(min_mn * 17.0 / 9.0) if max_mn >= mnthr1: if job == 'N': lwork = 130 * min_mn elif job == 'S': lwork = (min_mn + 130) * min_mn else: lwork = max( (min_mn + 130) * min_mn, (min_mn + 1) * min_mn + 32 * max_mn, ) else: lwork = 64 * (min_mn + max_mn) + 2 * min_mn else: mnthr = int(min_mn * 11.0 / 6.0) if m >= n: if m >= mnthr: if job == 'N': lwork = 131 * n elif job == 'S': lwork = max((131 + n) * n, (4 * n + 7) * n) else: lwork = max( (n + 131) * n, (n + 1) * n + 32 * m, (4 * n + 6) * n + m ) else: if job == 'N': lwork = 64 * m + 67 * n elif job == 'S': lwork = max(64 * m + 67 * n, (3 * n + 7) * n) else: lwork = (3 * n + 7) * n else: if n >= mnthr: if job == 'N': lwork = 131 * m elif job == 'S': lwork = max((m + 131) * m, (4 * m + 7) * m) else: lwork = max( (m + 131) * m, (m + 1) * m + 32 * n, (4 * m + 7) * m ) else: if job == 'N': lwork = 67 * m + 64 * n else: lwork = max(67 * m + 64 * n, (3 * m + 7) * m) s = nlcpy.empty((min_mn,) + a.shape[2:], dtype=f_dtype, order='F') work = nlcpy.empty(lwork, dtype=dtype) if a_complex: if job == 'N': lrwork = 7 * n * m else: lrwork = min_mn * max(5 * min_mn + 7, 2 * max(m, n) + 2 * min_mn + 1) else: lrwork = 1 rwork = nlcpy.empty(lrwork, dtype=f_dtype) iwork = nlcpy.empty(8 * min_mn, dtype='l') info = numpy.empty(1, dtype='l') fpe = request._get_fpe_flag() args = ( ord(job), a, s, u, vt, work, rwork, iwork, veo.OnStack(info, inout=veo.INTENT_OUT), veo.OnStack(fpe, inout=veo.INTENT_OUT), ) request._push_and_flush_request( 'nlcpy_svd', args, ) if c_order: s = nlcpy.asarray(nlcpy.moveaxis(s, 0, -1), order='C') else: s = nlcpy.moveaxis(s, 0, -1) if compute_uv: u = nlcpy.moveaxis(u, (1, 0), (-1, -2)) if not full_matrices: u = u[..., :m, :min_mn] if c_order: u = nlcpy.asarray(u, dtype=dtype, order='C') vt = nlcpy.asarray(nlcpy.moveaxis(vt, (1, 0), (-1, -2)), dtype, order='C') else: vt = nlcpy.moveaxis(nlcpy.asarray(vt, dtype), (1, 0), (-1, -2)) return u, s, vt else: return s
[ドキュメント]def cholesky(a): """Cholesky decomposition. Return the Cholesky decomposition, *L* * *L.H*, of the square matrix *a*, where *L* is lower-triangular and *.H* is the conjugate transpose operator (which is the ordinary transpose if *a* is real-valued). *a* must be Hermitian (symmetric if real-valued) and positive-definite. Only *L* is actually returned. Parameters ---------- a : (..., M, M) array_like Hermitian (symmetric if all elements are real), positive-definite input matrix. Returns ------- L : (..., M, M) ndarray Upper or lower-triangular Cholesky factor of *a*. Note ---- The Cholesky decomposition is often used as a fast way of solving :math:`Ax = b` (when *A* is both Hermitian/symmetric and positive-definite). First, we solve for y in :math:`Ly = b`, and then for x in :math:`L.Hx = y`. Examples -------- >>> import nlcpy as vp >>> A = vp.array([[1,-2j],[2j,5]]) >>> A array([[ 1.+0.j, -0.-2.j], [ 0.+2.j, 5.+0.j]]) >>> L = vp.linalg.cholesky(A) >>> L array([[1.+0.j, 0.+0.j], [0.+2.j, 1.+0.j]]) >>> vp.dot(L, vp.conjugate(L.T)) # verify that L * L.H = A array([[1.+0.j, 0.-2.j], [0.+2.j, 5.+0.j]]) """ a = nlcpy.asarray(a) util._assertRankAtLeast2(a) util._assertNdSquareness(a) if a.dtype == 'F': dtype = 'D' L_dtype = 'F' elif a.dtype == 'D': dtype = 'D' L_dtype = 'D' elif a.dtype == 'f': dtype = 'd' L_dtype = 'f' else: dtype = 'd' L_dtype = 'd' if a.size == 0: return nlcpy.empty(a.shape, dtype=L_dtype) # used to match the contiguous of result to numpy. c_order = a.flags.c_contiguous or sum([i > 1 for i in a.shape[:-2]]) < 2 a = nlcpy.array(nlcpy.moveaxis(a, (-1, -2), (1, 0)), dtype=dtype, order='F') info = numpy.empty(1, dtype='l') fpe = request._get_fpe_flag() args = ( a, veo.OnStack(info, inout=veo.INTENT_OUT), veo.OnStack(fpe, inout=veo.INTENT_OUT), ) request._push_and_flush_request( 'nlcpy_cholesky', args, callback=util._assertPositiveDefinite(info) ) if c_order: L = nlcpy.asarray(nlcpy.moveaxis(a, (1, 0), (-1, -2)), dtype=L_dtype, order='C') else: L = nlcpy.moveaxis(nlcpy.asarray(a, dtype=L_dtype), (1, 0), (-1, -2)) return L
[ドキュメント]def qr(a, mode='reduced'): """Computes the qr factorization of a matrix. Factor the matrix *a* as *qr*, where *q* is orthonormal and *r* is upper-triangular. Parameters ---------- a : (M, N) array_like Matrix to be factored. mode : {'reduced', 'complete', 'r', 'raw', 'full', 'economic'}, optional If K = min(M, N), then - 'reduced' : returns q, r with dimensions (M, K), (K, N) (default) - 'complete' : returns q, r with dimensions (M, M), (M, N) - 'r' : returns r only with dimensions (K, N) - 'raw' : returns h, tau with dimensions (N, M), (K,) - 'full' or 'f' : alias of 'reduced', deprecated - 'economic' or 'e' : returns h from 'raw', deprecated. Returns ------- q : ndarray, optional A matrix with orthonormal columns. When mode = 'complete' the result is an orthogonal/unitary matrix depending on whether or not a is real/complex. The determinant may be either +/- 1 in that case. r : ndarray, optional The upper-triangular matrix. (h, tau) : ndarray, optional The array h contains the Householder reflectors that generate q along with r. The tau array contains scaling factors for the reflectors. In the deprecated 'economic' mode only h is returned. Note ---- This is an interface to the LAPACK routines ``dgeqrf``, ``zgeqrf``, ``dorgqr``, and ``zungqr``. For more information on the qr factorization, see for example: https://en.wikipedia.org/wiki/QR_factorization Note that when 'raw' option is specified the returned arrays are of type "float64" or "complex128" and the h array is transposed to be FORTRAN compatible. Examples -------- >>> import numpy as np >>> import nlcpy as vp >>> from nlcpy import testing >>> a = vp.random.randn(9, 6) >>> q, r = vp.linalg.qr(a) >>> vp.testing.assert_allclose(a, vp.dot(q, r)) # a does equal qr >>> r2 = vp.linalg.qr(a, mode='r') >>> r3 = vp.linalg.qr(a, mode='economic') >>> # mode='r' returns the same r as mode='full' >>> vp.testing.assert_allclose(r, r2) >>> # But only triu parts are guaranteed equal when mode='economic' >>> vp.testing.assert_allclose(r, np.triu(r3[:6,:6], k=0)) Example illustrating a common use of qr: solving of least squares problems What are the least-squares-best *m* and *y0* in ``y = y0 + mx`` for the following data: {(0,1), (1,0), (1,2), (2,1)}. (Graph the points and you’ll see that it should be y0 = 0, m = 1.) The answer is provided by solving the over-determined matrix equation ``Ax = b``, where:: A = array([[0, 1], [1, 1], [1, 1], [2, 1]]) x = array([[y0], [m]]) b = array([[1], [0], [2], [1]]) If A = qr such that q is orthonormal (which is always possible via Gram-Schmidt), then ``x = inv(r) * (q.T) * b``. (In practice, however, we simply use :func:`lstsq`.) >>> A = vp.array([[0, 1], [1, 1], [1, 1], [2, 1]]) >>> A array([[0, 1], [1, 1], [1, 1], [2, 1]]) >>> b = vp.array([1, 0, 2, 1]) >>> q, r = vp.linalg.qr(A) >>> p = vp.dot(q.T, b) >>> vp.dot(vp.linalg.inv(r), p) array([1.1102230246251565e-16, 1.0000000000000002e+00]) """ if mode not in ('reduced', 'complete', 'r', 'raw'): if mode in ('f', 'full'): msg = "".join(( "The 'full' option is deprecated in favor of 'reduced'.\n", "For backward compatibility let mode default.")) warnings.warn(msg, DeprecationWarning, stacklevel=3) mode = 'reduced' elif mode in ('e', 'economic'): msg = "The 'economic' option is deprecated." warnings.warn(msg, DeprecationWarning, stacklevel=3) mode = 'economic' else: raise ValueError("Unrecognized mode '%s'" % mode) a = nlcpy.asarray(a) util._assertRank2(a) if a.dtype == 'F': dtype = 'D' a_dtype = 'F' elif a.dtype == 'D': dtype = 'D' a_dtype = 'D' elif a.dtype == 'f': dtype = 'd' a_dtype = 'f' else: dtype = 'd' a_dtype = 'd' m, n = a.shape if a.size == 0: if mode == 'reduced': return nlcpy.empty((m, 0), a_dtype), nlcpy.empty((0, n), a_dtype) elif mode == 'complete': return nlcpy.identity(m, a_dtype), nlcpy.empty((m, n), a_dtype) elif mode == 'r': return nlcpy.empty((0, n), a_dtype) elif mode == 'raw': return nlcpy.empty((n, m), dtype), nlcpy.empty((0,), dtype) else: return nlcpy.empty((m, n), a_dtype), nlcpy.empty((0,), a_dtype) a = nlcpy.asarray(a, dtype=dtype, order='F') k = min(m, n) if mode == 'complete': if m > n: x = nlcpy.empty((m, m), dtype=dtype, order='F') x[:m, :n] = a a = x r_shape = (m, n) elif mode in ('r', 'reduced', 'economic'): r_shape = (k, n) else: r_shape = 1 jobq = 0 if mode in ('r', 'raw', 'economic') else 1 tau = nlcpy.empty(k, dtype=dtype) r = nlcpy.zeros(r_shape, dtype=dtype) work = nlcpy.empty(n * 64, dtype=dtype) fpe = request._get_fpe_flag() args = ( m, n, jobq, a, tau, r, work, veo.OnStack(fpe, inout=veo.INTENT_OUT), ) request._push_and_flush_request( 'nlcpy_qr', args, ) if mode == 'raw': return a.T, tau if mode == 'r': return nlcpy.asarray(r, dtype=a_dtype) if mode == 'economic': return nlcpy.asarray(a, dtype=a_dtype) mc = m if mode == 'complete' else k q = nlcpy.asarray(a[:, :mc], dtype=a_dtype, order='C') r = nlcpy.asarray(r, dtype=a_dtype, order='C') return q, r