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

#
# * The source code in this file is developed independently by NEC Corporation.
#
# # 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.
#

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


[ドキュメント]def solve(a, b): """Solves a linear matrix equation, or system of linear scalar equations. Computes the "exact" solution, *x*, of the well-determined, i.e., full rank, linear matrix equation :math:`ax = b`. Parameters ---------- a : (..., M, M) array_like Coefficient matrix. b : {(..., M,), (..., M, K)} array_like Ordinate or "dependent variable" values. Returns ------- x : {(..., M,), (..., M, K)} ndarray Solution to the system a x = b. Returned shape is identical to *b*. Note ---- The solutions are computed using LAPACK routine ``_gesv``. `a` must be square and of full-rank, i.e., all rows (or, equivalently, columns) must be linearly independent; if either is not true, use :func:`lstsq` for the least-squares best "solution" of the system/equation. Examples -------- Solve the system of equations ``3 * x0 + x1 = 9`` and ``x0 + 2 * x1 = 8``: >>> import nlcpy as vp >>> a = vp.array([[3,1], [1,2]]) >>> b = vp.array([9,8]) >>> x = vp.linalg.solve(a, b) >>> x array([2., 3.]) """ a = nlcpy.asarray(a) b = nlcpy.asarray(b) c_order = (a.flags.c_contiguous or a.ndim < 4 or a.ndim - b.ndim < 2 and b.flags.c_contiguous) and \ not (a.ndim < b.ndim and not b.flags.c_contiguous) util._assertRankAtLeast2(a) util._assertNdSquareness(a) if a.ndim - 1 == b.ndim: if a.shape[-1] != b.shape[-1]: raise ValueError( 'solve1: Input operand 1 has a mismatch in ' 'its core dimension 0, with gufunc signature (m,m),(m)->(m) ' '(size {0} is different from {1})'.format(b.shape[-1], a.shape[-1])) elif b.ndim == 1: raise ValueError( 'solve: Input operand 1 does not have enough dimensions ' '(has 1, gufunc core with signature (m,m),(m,n)->(m,n) requires 2)') else: if a.shape[-1] != b.shape[-2]: raise ValueError( 'solve: Input operand 1 has a mismatch in ' 'its core dimension 0, with gufunc signature (m,m),(m,n)->(m,n) ' '(size {0} is different from {1})'.format(b.shape[-2], a.shape[-1])) if b.ndim == 1 or a.ndim - 1 == b.ndim and a.shape[-1] == b.shape[-1]: tmp = 1 _newaxis = (None, ) else: tmp = 2 _newaxis = (None, None) for i in range(1, min(a.ndim - 2, b.ndim - tmp) + 1): if a.shape[-2 - i] != b.shape[-tmp - i] and \ 1 not in (a.shape[-2 - i], b.shape[-tmp - i]): raise ValueError( 'operands could not be broadcast together with ' 'remapped shapes [original->remapped]: {0}->({1}) ' '{2}->({3}) and requested shape ({4})'.format( str(a.shape).replace(' ', ''), str(a.shape[:-2] + _newaxis).replace(' ', ''). replace('None', 'newaxis').strip('(,)'), str(b.shape).replace(' ', ''), str(b.shape[:-tmp] + _newaxis).replace(' ', ''). replace('None', 'newaxis'). replace('None', 'newaxis').strip('(,)'), str(b.shape[-tmp:]).replace(' ', '').strip('(,)') )) if a.dtype.char in 'FD' or b.dtype.char in 'FD': dtype = 'complex128' if a.dtype.char in 'fF' and b.dtype.char in 'fF': x_dtype = 'complex64' else: x_dtype = 'complex128' else: dtype = 'float64' if a.dtype.char == 'f' and b.dtype.char == 'f': x_dtype = 'float32' else: x_dtype = 'float64' x_shape = b.shape if b.ndim == a.ndim - 1: b = b[..., nlcpy.newaxis] diff = abs(a.ndim - b.ndim) if a.ndim < b.ndim: bcast_shape = [b.shape[i] if b.shape[i] != 1 or i < diff else a.shape[i - diff] for i in range(b.ndim - 2)] else: bcast_shape = [a.shape[i] if a.shape[i] != 1 or i < diff else b.shape[i - diff] for i in range(a.ndim - 2)] bcast_shape_a = bcast_shape + list(a.shape[-2:]) bcast_shape_b = bcast_shape + list(b.shape[-2:]) a = nlcpy.broadcast_to(a, bcast_shape_a) if bcast_shape_b != list(b.shape): b = nlcpy.broadcast_to(b, bcast_shape_b) x_shape = b.shape if b.size == 0: return nlcpy.empty(x_shape, dtype=x_dtype) a = nlcpy.array(nlcpy.moveaxis(a, (-1, -2), (1, 0)), dtype=dtype, order='F') b = nlcpy.array(nlcpy.moveaxis(b, (-1, -2), (1, 0)), dtype=dtype, order='F') info = numpy.empty(1, dtype='l') fpe = request._get_fpe_flag() args = ( a, b, veo.OnStack(info, inout=veo.INTENT_OUT), veo.OnStack(fpe, inout=veo.INTENT_OUT), ) request._push_and_flush_request( 'nlcpy_solve', args, callback=util._assertNotSingular(info) ) if c_order: x = nlcpy.moveaxis(b, (1, 0), (-1, -2)).reshape(x_shape) return nlcpy.asarray(x, x_dtype, 'C') else: x = nlcpy.asarray(b, x_dtype) return nlcpy.moveaxis(x, (1, 0), (-1, -2)).reshape(x_shape)
def _lstsq_errchk(info): def _info_check(*args): if info > 0: raise util.LinAlgError('a singular value did not converge') return _info_check
[ドキュメント]def lstsq(a, b, rcond='warn'): """Returns the least-squares solution to a linear matrix equation. Solves the equation :math:`ax = b` by computing a vector *x* that minimizes the squared Euclidean 2-norm :math:`|b-ax|^2_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. Parameters ---------- a : (M,N) array_like "Coefficient" matrix. b : {(M,), (M, K)} array_like Ordinate or "dependent variable" values. If *b* is two-dimensional, the least-squares solution is calculated for each of the *K* columns of *b*. rcond : float, optional Cut-off ratio for small singular values of *a*. For the purposes of rank determination, singular values are treated as zero if they are smaller than *rcond* times the largest singular value of *a*. Returns ------- x : {(N,), (N, K)} ndarray Least-squares solution. If *b* is two-dimensional, the solutions are in the *K* columns of *x*. residuals : {(1,), (K,), (0,)} ndarray Sums of residuals; squared Euclidean 2-norm for each column in ``b - a@x``. If the rank of *a* is < N or M <= N, this is an empty array. If *b* is 1-dimensional, this is a (1,) shape array. Otherwise the shape is (K,). rank : *int* Rank of matrix *a*. s : (min(M, N),) ndarray Singular values of *a*. Note ---- If `b` is a matrix, then all array results are returned as matrices. Examples -------- .. plot:: :align: center Fit a line, ``y = mx + c``, through some noisy data-points: >>> import nlcpy as vp >>> x = vp.array([0, 1, 2, 3]) >>> y = vp.array([-1, 0.2, 0.9, 2.1]) By examining the coefficients, we see that the line should have a gradient of roughly 1 and cut the y-axis at, more or less, -1. We can rewrite the line equation as ``y = Ap``, where ``A = [[x 1]]`` and ``p = [[m], [c]]``. Now use lstsq to solve for *p*: >>> A = vp.array((x, vp.ones(len(x)))).T >>> A array([[0., 1.], [1., 1.], [2., 1.], [3., 1.]]) >>> m, c = vp.linalg.lstsq(A, y, rcond=None)[0] >>> m, c # doctest: +SKIP (array(1.), array(-0.95)) # may vary Plot the data along with the fitted line: >>> import matplotlib.pyplot as plt >>> _ = plt.plot(x, y, 'o', label='Original data', markersize=15) >>> _ = plt.plot(x, m*x + c, 'r', label='Fitted line') >>> _ = plt.legend() >>> plt.show() """ a = nlcpy.asarray(a) b = nlcpy.asarray(b) b_ndim = b.ndim if b_ndim == 1: b = b[:, nlcpy.newaxis] util._assertRank2(a, b) if a.shape[-2] != b.shape[-2]: raise util.LinAlgError('Incompatible dimensions') a_complex = a.dtype.char in 'FD' b_complex = b.dtype.char in 'FD' if a_complex or b_complex: if a.dtype.char in 'fF' and b.dtype.char in 'fF': x_dtype = 'F' f_dtype = 'f' else: x_dtype = 'D' f_dtype = 'd' else: if a.dtype.char == 'f' and b.dtype.char == 'f': x_dtype = 'f' f_dtype = 'f' else: x_dtype = 'd' f_dtype = 'd' m = a.shape[-2] n = a.shape[-1] k = b.shape[-1] minmn = min(m, n) maxmn = max(m, n) k_extend = (k == 0) if k_extend: b = nlcpy.zeros([m, 1], dtype=b.dtype) k = 1 if minmn == 0: if n > 0: x = nlcpy.zeros([n, k], dtype=x_dtype) residuals = nlcpy.array([], dtype=f_dtype) else: x = nlcpy.array([], dtype=x_dtype) if b_complex: br = nlcpy.asarray(b.real, dtype='d') bi = nlcpy.asarray(b.imag, dtype='d') square_b = nlcpy.multiply(br, br) + nlcpy.multiply(bi, bi) else: square_b = nlcpy.multiply(b, b, dtype='d') residuals = nlcpy.add.reduce(square_b, dtype=f_dtype) return (x, residuals, 0, nlcpy.array([], dtype=f_dtype)) if rcond == 'warn': rcond = -1.0 if rcond is None: rcond = nlcpy.finfo(x_dtype).eps * max(m, n) rcond = nlcpy.asarray(rcond, dtype=f_dtype) a = nlcpy.array(a, dtype=x_dtype, order='F') if m < n: _b = nlcpy.empty([n, k], dtype=x_dtype, order='F') _b[:m, :] = b b = _b else: b = nlcpy.array(b, dtype=x_dtype, order='F') s = nlcpy.empty(min(m, n), dtype=f_dtype) rank = numpy.empty(1, dtype='l') nlvl = max(0, int(nlcpy.log(minmn / 26.0) / nlcpy.log(2)) + 1) mnthr = int(minmn * 1.6) mm = m lwork = 1 if a_complex or b_complex: _tmp = 2 wlalsd = minmn * k lrwork = 10 * maxmn + 2 * maxmn * 25 + 8 * maxmn * nlvl + 3 * 25 * k \ + max(26 ** 2, n * (1 + k) + 2 * k) else: _tmp = 3 wlalsd = 9 * minmn + 2 * minmn * 25 + 8 * minmn * nlvl + minmn * k + 26 ** 2 lrwork = 1 if m >= n: if m >= mnthr: mm = n lwork = max(65 * n, n + 32 * k) lwork = max(lwork, _tmp * n + (mm + n) * 64, _tmp * n + 32 * k, _tmp * n + 32 * n - 32, _tmp * n + wlalsd) else: if n >= mnthr: lwork = max(lwork, m * m + 4 * m + wlalsd, m * m + 132 * m, m * m + 4 * m + 32 * k, m * m + (k + 1) * m, m * m + n + m) else: lwork = max(lwork, 3 * m + wlalsd, 3 * m + 64 * (n + m), 3 * m + 32 * k) liwork = minmn * (11 + 3 * nlvl) work = nlcpy.empty(lwork, dtype=f_dtype) iwork = nlcpy.empty(liwork, dtype='l') rwork = nlcpy.empty(lrwork, dtype=f_dtype) info = numpy.empty(1, dtype='l') fpe = request._get_fpe_flag() args = ( a, b, s, work, rwork, iwork, rcond, veo.OnStack(rank, inout=veo.INTENT_OUT), veo.OnStack(info, inout=veo.INTENT_OUT), veo.OnStack(fpe, inout=veo.INTENT_OUT), ) request._push_and_flush_request( 'nlcpy_lstsq', args, callback=_lstsq_errchk(info), sync=True ) if rank < n or m <= n: residuals = nlcpy.array([], dtype=f_dtype) else: _b = b[n:] square_b = nlcpy.multiply(_b.real, _b.real) + nlcpy.multiply(_b.imag, _b.imag) residuals = nlcpy.add.reduce(square_b, dtype=f_dtype) if k_extend: x = b[..., :0] residuals = residuals[..., :0] elif b_ndim == 1: x = nlcpy.asarray(b[:n, 0], dtype=x_dtype, order='C') else: x = nlcpy.asarray(b[:n, :], dtype=x_dtype, order='C') return (x, residuals, rank[0], s)
[ドキュメント]def inv(a): """Computes the (multiplicative) inverse of a matrix. Given a square matrix *a*, return the matrix *ainv* satisfying :: dot(a, ainv) = dot(ainv, a) = eye(a.shape[0]). Parameters ---------- a : (..., M, M) array_like Matrix to be inverted. Returns ------- ainv : (..., M, M) ndarray (Multiplicative) inverse of the matrix *a*. Note ---- Broadcasting rules apply, see the :ref:`nlcpy.linalg <nlcpy_linalg>` documentation for details. Examples -------- >>> import nlcpy as vp >>> from nlcpy import testing >>> a = vp.array([[1., 2.], [3., 4.]]) >>> ainv = vp.linalg.inv(a) >>> vp.testing.assert_allclose(vp.dot(a, ainv), vp.eye(2), atol=1e-8, rtol=1e-5) >>> vp.testing.assert_allclose(vp.dot(ainv, a), vp.eye(2), atol=1e-8, rtol=1e-5) Inverses of several matrices can be computed at once: >>> a = vp.array([[[1., 2.], [3., 4.]], [[1, 3], [3, 5]]]) >>> vp.linalg.inv(a) # doctest: +SKIP array([[[-2. , 1. ], [ 1.5 , -0.5 ]], <BLANKLINE> [[-1.25, 0.75], [ 0.75, -0.25]]]) """ a = nlcpy.asarray(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 util._assertRankAtLeast2(a) util._assertNdSquareness(a) if a.dtype.char in 'FD': dtype = 'D' if a.dtype.char in 'fF': ainv_dtype = 'F' else: ainv_dtype = 'D' else: dtype = 'd' if a.dtype.char == 'f': ainv_dtype = 'f' else: ainv_dtype = 'd' if a.size == 0: return nlcpy.asarray(a, dtype=ainv_dtype) a = nlcpy.array(nlcpy.moveaxis(a, (-1, -2), (1, 0)), dtype=dtype, order='F') ipiv = nlcpy.empty(a.shape[-1], dtype='l') work = nlcpy.empty(a.shape[-1] * 256) info = numpy.empty(1, dtype='l') fpe = request._get_fpe_flag() args = ( a, ipiv, work, veo.OnStack(info, inout=veo.INTENT_OUT), veo.OnStack(fpe, inout=veo.INTENT_OUT), ) request._push_and_flush_request( 'nlcpy_inv', args, callback=util._assertNotSingular(info) ) if c_order: a = nlcpy.moveaxis(a, (1, 0), (-1, -2)) return nlcpy.asarray(a, dtype=ainv_dtype, order='C') else: a = nlcpy.asarray(a, dtype=ainv_dtype) return nlcpy.moveaxis(a, (1, 0), (-1, -2))