"""QR decomposition functions.""" import numpy as np from scipy._lib._util import _apply_over_batch, _deprecate_dtypes from scipy._lib.deprecation import _NoValue import warnings # Local imports from .lapack import _normalize_lapack_dtype, get_lapack_funcs, HAS_ILP64 from ._misc import _datacopied from ._basic import _format_emit_errors_warnings from . import _batched_linalg __all__ = ['qr', 'qr_multiply', 'rq'] def safecall(f, name, *args, **kwargs): """Call a LAPACK routine, determining lwork automatically and handling error return values""" lwork = kwargs.get("lwork", None) if lwork in (None, -1): kwargs['lwork'] = -1 ret = f(*args, **kwargs) kwargs['lwork'] = ret[-2][0].real.astype(np.int_) ret = f(*args, **kwargs) if ret[-1] < 0: raise ValueError(f"illegal value in {-ret[-1]}th argument of internal {name}") return ret[:-2] def qr(a, overwrite_a=False, lwork=_NoValue, mode="full", pivoting=False, check_finite=True): """ Compute QR decomposition of a matrix. Calculate the decomposition ``A = Q R`` where Q is unitary/orthogonal and R upper triangular. Parameters ---------- a : (..., M, N) array_like Matrix to be decomposed overwrite_a : bool, optional Whether to overwrite data in `a` (may improve performance). Default is False. See :ref:`tutorial_linalg_overwrite` for details. lwork : int, optional Work array size, lwork >= a.shape[1]. If None or -1, an optimal size is computed. .. deprecated:: 1.18.0 This keyword is deprecated as well as no longer in use and will be removed in 1.20.0. mode : {'full', 'r', 'economic', 'raw'}, optional Determines what information is to be returned: either both Q and R ('full', default), only R ('r') or both Q and R but computed in economy-size ('economic', see Notes). The final option 'raw' (added in SciPy 0.11) makes the function return two matrices (Q, TAU) in the internal format used by LAPACK. pivoting : bool, optional Whether or not factorization should include pivoting for rank-revealing qr decomposition. If pivoting, compute the decomposition ``A[..., :, P] = Q @ R`` as above, but where P is chosen such that the diagonal of R is non-increasing. Equivalently, albeit less efficiently, an explicit P matrix may be formed explicitly by permuting the rows or columns (depending on the side of the equation on which it is to be used) of an identity matrix. See Examples. check_finite : bool, optional Whether to check that the input matrix contains only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs. Returns ------- Q : float or complex ndarray Of shape (..., M, M), or (..., M, K) for ``mode='economic'``. Not returned if ``mode='r'``. Replaced by tuple ``(Q, TAU)`` if ``mode='raw'``. R : float or complex ndarray Of shape (..., M, N), or (..., K, N) for ``mode in ['economic', 'raw']``. ``K = min(M, N)``. P : int ndarray Of shape (..., N,) for ``pivoting=True``. Not returned if ``pivoting=False``. Raises ------ LinAlgError Raised if decomposition fails Notes ----- This is an interface to the LAPACK routines dgeqrf, zgeqrf, dorgqr, zungqr, dgeqp3, and zgeqp3. If ``mode=economic``, the shapes of Q and R are (..., M, K) and (..., K, N) instead of (..., M,M) and (..., M,N), with ``K=min(M,N)``. Examples -------- >>> import numpy as np >>> from scipy import linalg >>> rng = np.random.default_rng() >>> a = rng.standard_normal((9, 6)) >>> q, r = linalg.qr(a) >>> np.allclose(a, np.dot(q, r)) True >>> q.shape, r.shape ((9, 9), (9, 6)) >>> r2 = linalg.qr(a, mode='r') >>> np.allclose(r, r2) True >>> q3, r3 = linalg.qr(a, mode='economic') >>> q3.shape, r3.shape ((9, 6), (6, 6)) >>> q4, r4, p4 = linalg.qr(a, pivoting=True) >>> d = np.abs(np.diag(r4)) >>> np.all(d[1:] <= d[:-1]) True >>> np.allclose(a[:, p4], np.dot(q4, r4)) True >>> P = np.eye(p4.size)[p4] >>> np.allclose(a, np.dot(q4, r4) @ P) True >>> np.allclose(a @ P.T, np.dot(q4, r4)) True >>> q4.shape, r4.shape, p4.shape ((9, 9), (9, 6), (6,)) >>> q5, r5, p5 = linalg.qr(a, mode='economic', pivoting=True) >>> q5.shape, r5.shape, p5.shape ((9, 6), (6, 6), (6,)) >>> P = np.eye(6)[:, p5] >>> np.allclose(a @ P, np.dot(q5, r5)) True """ # structure mappings, keep in sync with the C side modes = { "full": 1, "qr": 1, # equivalent to `full` "r": 11, "raw": 21, "economic": 31 } # 'qr' was the old default, equivalent to 'full'. Neither 'full' nor # 'qr' are used below. # 'raw' is used internally by qr_multiply if mode not in modes.keys(): raise ValueError(f"Mode argument should be one of {list(modes.keys())}") modeFlag = modes[mode] # convert the string to an int for the C enum if check_finite: a1 = np.asarray_chkfinite(a) else: a1 = np.asarray(a) _deprecate_dtypes("linalg.qr", a1) if a1.ndim < 2: raise ValueError("Expected at least a 2-D array") M, N = a1.shape[-2], a1.shape[-1] # Throw error for backwards compat, else raise DeprecationWarning if lwork is not _NoValue: if lwork is not None and lwork != -1 and lwork <= M: raise ValueError(f"lwork should be None, -1 or > M, got {lwork}") else: warnings.warn( "scipy.linalg: the `lwork` keyword is deprecated and no longer in use" " as of SciPy 1.18.0 and will be removed in SciPy 1.20.0", DeprecationWarning, stacklevel=2 ) # First normalize dtypes to ensure consistent return types a1, overwrite_a = _normalize_lapack_dtype(a1, overwrite_a) # accommodate empty arrays if a1.size == 0: K = min(M, N) batch_shape = a1.shape[:-2] if mode not in ['economic', 'raw']: Q = np.empty_like(a1, shape=batch_shape + (M, M)) Q[..., :, :] = np.identity(M) R = np.empty_like(a1) else: Q = np.empty_like(a1, shape=batch_shape + (M, K)) R = np.empty_like(a1, shape=batch_shape + (K, N)) if pivoting: Rj = R, np.arange(N, dtype=np.int64 if HAS_ILP64 else np.int32) else: Rj = R, if mode == 'r': return Rj elif mode == 'raw': qr = np.empty_like(a1, shape=batch_shape + (M, N)) tau = np.zeros_like(a1, shape=batch_shape + (K,)) return ((qr, tau),) + Rj return (Q,) + Rj if not (a1.flags['ALIGNED'] and a1.dtype.byteorder == '='): overwrite_a = True a1 = a1.copy() overwrite_a = overwrite_a or (_datacopied(a1, a)) overwrite_a = overwrite_a and (a1.ndim == 2) and a1.flags["F_CONTIGUOUS"] # heavy lifting Q, R, tau, jpvt, err_lst = _batched_linalg._qr(a1, overwrite_a, modeFlag, pivoting) if err_lst: _format_emit_errors_warnings(err_lst) # construct return objects if pivoting: Rj = R, jpvt else: Rj = (R,) if modeFlag == modes["raw"]: Q = (Q, tau) elif modeFlag == modes["economic"] and M < N: Q = Q[..., :, :M] if modeFlag == modes["r"]: return Rj else: return (Q,) + Rj @_apply_over_batch(('a', 2), ('c', '1|2')) def qr_multiply(a, c, mode='right', pivoting=False, conjugate=False, overwrite_a=False, overwrite_c=False): """ Calculate the QR decomposition and multiply Q with a matrix. Calculate the decomposition ``A = Q R`` where Q is unitary/orthogonal and R upper triangular. Multiply Q with a vector or a matrix c. Parameters ---------- a : (M, N), array_like Input array c : array_like Input array to be multiplied by ``q``. mode : {'left', 'right'}, optional ``Q @ c`` is returned if mode is 'left', ``c @ Q`` is returned if mode is 'right'. The shape of c must be appropriate for the matrix multiplications, if mode is 'left', ``min(a.shape) == c.shape[0]``, if mode is 'right', ``a.shape[0] == c.shape[1]``. pivoting : bool, optional Whether or not factorization should include pivoting for rank-revealing qr decomposition, see the documentation of qr. conjugate : bool, optional Whether Q should be complex-conjugated. This might be faster than explicit conjugation. overwrite_a : bool, optional Whether to overwrite data in `a` (may improve performance). Default is False. See :ref:`tutorial_linalg_overwrite` for details. overwrite_c : bool, optional Whether data in c is overwritten (may improve performance). If this is used, c must be big enough to keep the result, i.e. ``c.shape[0]`` = ``a.shape[0]`` if mode is 'left'. See :ref:`tutorial_linalg_overwrite` for details. Returns ------- CQ : ndarray The product of ``Q`` and ``c``. R : (K, N), ndarray R array of the resulting QR factorization where ``K = min(M, N)``. P : (N,) ndarray Integer pivot array. Only returned when ``pivoting=True``. Raises ------ LinAlgError Raised if QR decomposition fails. Notes ----- This is an interface to the LAPACK routines ``?GEQRF``, ``?ORMQR``, ``?UNMQR``, and ``?GEQP3``. .. versionadded:: 0.11.0 Examples -------- >>> import numpy as np >>> from scipy.linalg import qr_multiply, qr >>> A = np.array([[1, 3, 3], [2, 3, 2], [2, 3, 3], [1, 3, 2]]) >>> qc, r1, piv1 = qr_multiply(A, 2*np.eye(4), pivoting=1) >>> qc array([[-1., 1., -1.], [-1., -1., 1.], [-1., -1., -1.], [-1., 1., 1.]]) >>> r1 array([[-6., -3., -5. ], [ 0., -1., -1.11022302e-16], [ 0., 0., -1. ]]) >>> piv1 array([1, 0, 2], dtype=int32) >>> q2, r2, piv2 = qr(A, mode='economic', pivoting=1) >>> np.allclose(2*q2 - qc, np.zeros((4, 3))) True """ if mode not in ['left', 'right']: raise ValueError("Mode argument can only be 'left' or 'right' but " f"not '{mode}'") c = np.asarray_chkfinite(c) if c.ndim < 2: onedim = True c = np.atleast_2d(c) if mode == "left": c = c.T else: onedim = False a = np.atleast_2d(np.asarray(a)) # chkfinite done in qr M, N = a.shape if mode == 'left': if c.shape[0] != min(M, N + overwrite_c*(M-N)): raise ValueError('Array shapes are not compatible for Q @ c' f' operation: {a.shape} vs {c.shape}') else: if M != c.shape[1]: raise ValueError('Array shapes are not compatible for c @ Q' f' operation: {c.shape} vs {a.shape}') raw = qr(a, overwrite_a, mode="raw", pivoting=pivoting) Q, tau = raw[0] # accommodate empty arrays if c.size == 0: return (np.empty_like(c),) + raw[1:] gor_un_mqr, = get_lapack_funcs(('ormqr',), (Q,)) if gor_un_mqr.typecode in ('s', 'd'): trans = "T" else: trans = "C" Q = Q[:, :min(M, N)] if M > N and mode == "left" and not overwrite_c: if conjugate: cc = np.zeros((c.shape[1], M), dtype=c.dtype, order="F") cc[:, :N] = c.T else: cc = np.zeros((M, c.shape[1]), dtype=c.dtype, order="F") cc[:N, :] = c trans = "N" if conjugate: lr = "R" else: lr = "L" overwrite_c = True elif c.flags["C_CONTIGUOUS"] and trans == "T" or conjugate: cc = c.T if mode == "left": lr = "R" else: lr = "L" else: trans = "N" cc = c if mode == "left": lr = "L" else: lr = "R" cQ, = safecall(gor_un_mqr, "gormqr/gunmqr", lr, trans, Q, tau, cc, overwrite_c=overwrite_c) if trans != "N": cQ = cQ.T if mode == "right": cQ = cQ[:, :min(M, N)] if onedim: cQ = cQ.ravel() return (cQ,) + raw[1:] @_apply_over_batch(('a', 2)) def rq(a, overwrite_a=False, lwork=None, mode='full', check_finite=True): """ Compute RQ decomposition of a matrix. Calculate the decomposition ``A = R Q`` where Q is unitary/orthogonal and R upper triangular. Parameters ---------- a : (M, N) array_like Matrix to be decomposed overwrite_a : bool, optional Whether to overwrite data in `a` (may improve performance). Default is False. See :ref:`tutorial_linalg_overwrite` for details. lwork : int, optional Work array size, lwork >= a.shape[1]. If None or -1, an optimal size is computed. mode : {'full', 'r', 'economic'}, optional Determines what information is to be returned: either both Q and R ('full', default), only R ('r') or both Q and R but computed in economy-size ('economic', see Notes). check_finite : bool, optional Whether to check that the input matrix contains only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs. Returns ------- R : float or complex ndarray Of shape (M, N) or (M, K) for ``mode='economic'``. ``K = min(M, N)``. Q : float or complex ndarray Of shape (N, N) or (K, N) for ``mode='economic'``. Not returned if ``mode='r'``. Raises ------ LinAlgError If decomposition fails. Notes ----- This is an interface to the LAPACK routines sgerqf, dgerqf, cgerqf, zgerqf, sorgrq, dorgrq, cungrq and zungrq. If ``mode=economic``, the shapes of Q and R are (K, N) and (M, K) instead of (N,N) and (M,N), with ``K=min(M,N)``. Examples -------- >>> import numpy as np >>> from scipy import linalg >>> rng = np.random.default_rng() >>> a = rng.standard_normal((6, 9)) >>> r, q = linalg.rq(a) >>> np.allclose(a, r @ q) True >>> r.shape, q.shape ((6, 9), (9, 9)) >>> r2 = linalg.rq(a, mode='r') >>> np.allclose(r, r2) True >>> r3, q3 = linalg.rq(a, mode='economic') >>> r3.shape, q3.shape ((6, 6), (6, 9)) """ if mode not in ['full', 'r', 'economic']: raise ValueError( "Mode argument should be one of ['full', 'r', 'economic']") if check_finite: a1 = np.asarray_chkfinite(a) else: a1 = np.asarray(a) if len(a1.shape) != 2: raise ValueError('expected matrix') M, N = a1.shape # accommodate empty arrays if a1.size == 0: K = min(M, N) if not mode == 'economic': R = np.empty_like(a1) Q = np.empty_like(a1, shape=(N, N)) Q[...] = np.identity(N) else: R = np.empty_like(a1, shape=(M, K)) Q = np.empty_like(a1, shape=(K, N)) if mode == 'r': return R return R, Q overwrite_a = overwrite_a or (_datacopied(a1, a)) gerqf, = get_lapack_funcs(('gerqf',), (a1,)) rq, tau = safecall(gerqf, 'gerqf', a1, lwork=lwork, overwrite_a=overwrite_a) if not mode == 'economic' or N < M: R = np.triu(rq, N-M) else: R = np.triu(rq[-M:, -M:]) if mode == 'r': return R gor_un_grq, = get_lapack_funcs(('orgrq',), (rq,)) if N < M: Q, = safecall(gor_un_grq, "gorgrq/gungrq", rq[-N:], tau, lwork=lwork, overwrite_a=1) elif mode == 'economic': Q, = safecall(gor_un_grq, "gorgrq/gungrq", rq, tau, lwork=lwork, overwrite_a=1) else: rq1 = np.empty((N, N), dtype=rq.dtype) rq1[-M:] = rq Q, = safecall(gor_un_grq, "gorgrq/gungrq", rq1, tau, lwork=lwork, overwrite_a=1) return R, Q