Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/scipy/linalg/_decomp_qr.py: 6%
129 statements
« prev ^ index » next coverage.py v7.3.1, created at 2023-09-23 06:43 +0000
1"""QR decomposition functions."""
2import numpy
4# Local imports
5from .lapack import get_lapack_funcs
6from ._misc import _datacopied
8__all__ = ['qr', 'qr_multiply', 'rq']
11def safecall(f, name, *args, **kwargs):
12 """Call a LAPACK routine, determining lwork automatically and handling
13 error return values"""
14 lwork = kwargs.get("lwork", None)
15 if lwork in (None, -1):
16 kwargs['lwork'] = -1
17 ret = f(*args, **kwargs)
18 kwargs['lwork'] = ret[-2][0].real.astype(numpy.int_)
19 ret = f(*args, **kwargs)
20 if ret[-1] < 0:
21 raise ValueError("illegal value in %dth argument of internal %s"
22 % (-ret[-1], name))
23 return ret[:-2]
26def qr(a, overwrite_a=False, lwork=None, mode='full', pivoting=False,
27 check_finite=True):
28 """
29 Compute QR decomposition of a matrix.
31 Calculate the decomposition ``A = Q R`` where Q is unitary/orthogonal
32 and R upper triangular.
34 Parameters
35 ----------
36 a : (M, N) array_like
37 Matrix to be decomposed
38 overwrite_a : bool, optional
39 Whether data in `a` is overwritten (may improve performance if
40 `overwrite_a` is set to True by reusing the existing input data
41 structure rather than creating a new one.)
42 lwork : int, optional
43 Work array size, lwork >= a.shape[1]. If None or -1, an optimal size
44 is computed.
45 mode : {'full', 'r', 'economic', 'raw'}, optional
46 Determines what information is to be returned: either both Q and R
47 ('full', default), only R ('r') or both Q and R but computed in
48 economy-size ('economic', see Notes). The final option 'raw'
49 (added in SciPy 0.11) makes the function return two matrices
50 (Q, TAU) in the internal format used by LAPACK.
51 pivoting : bool, optional
52 Whether or not factorization should include pivoting for rank-revealing
53 qr decomposition. If pivoting, compute the decomposition
54 ``A P = Q R`` as above, but where P is chosen such that the diagonal
55 of R is non-increasing.
56 check_finite : bool, optional
57 Whether to check that the input matrix contains only finite numbers.
58 Disabling may give a performance gain, but may result in problems
59 (crashes, non-termination) if the inputs do contain infinities or NaNs.
61 Returns
62 -------
63 Q : float or complex ndarray
64 Of shape (M, M), or (M, K) for ``mode='economic'``. Not returned
65 if ``mode='r'``. Replaced by tuple ``(Q, TAU)`` if ``mode='raw'``.
66 R : float or complex ndarray
67 Of shape (M, N), or (K, N) for ``mode in ['economic', 'raw']``.
68 ``K = min(M, N)``.
69 P : int ndarray
70 Of shape (N,) for ``pivoting=True``. Not returned if
71 ``pivoting=False``.
73 Raises
74 ------
75 LinAlgError
76 Raised if decomposition fails
78 Notes
79 -----
80 This is an interface to the LAPACK routines dgeqrf, zgeqrf,
81 dorgqr, zungqr, dgeqp3, and zgeqp3.
83 If ``mode=economic``, the shapes of Q and R are (M, K) and (K, N) instead
84 of (M,M) and (M,N), with ``K=min(M,N)``.
86 Examples
87 --------
88 >>> import numpy as np
89 >>> from scipy import linalg
90 >>> rng = np.random.default_rng()
91 >>> a = rng.standard_normal((9, 6))
93 >>> q, r = linalg.qr(a)
94 >>> np.allclose(a, np.dot(q, r))
95 True
96 >>> q.shape, r.shape
97 ((9, 9), (9, 6))
99 >>> r2 = linalg.qr(a, mode='r')
100 >>> np.allclose(r, r2)
101 True
103 >>> q3, r3 = linalg.qr(a, mode='economic')
104 >>> q3.shape, r3.shape
105 ((9, 6), (6, 6))
107 >>> q4, r4, p4 = linalg.qr(a, pivoting=True)
108 >>> d = np.abs(np.diag(r4))
109 >>> np.all(d[1:] <= d[:-1])
110 True
111 >>> np.allclose(a[:, p4], np.dot(q4, r4))
112 True
113 >>> q4.shape, r4.shape, p4.shape
114 ((9, 9), (9, 6), (6,))
116 >>> q5, r5, p5 = linalg.qr(a, mode='economic', pivoting=True)
117 >>> q5.shape, r5.shape, p5.shape
118 ((9, 6), (6, 6), (6,))
120 """
121 # 'qr' was the old default, equivalent to 'full'. Neither 'full' nor
122 # 'qr' are used below.
123 # 'raw' is used internally by qr_multiply
124 if mode not in ['full', 'qr', 'r', 'economic', 'raw']:
125 raise ValueError("Mode argument should be one of ['full', 'r',"
126 "'economic', 'raw']")
128 if check_finite:
129 a1 = numpy.asarray_chkfinite(a)
130 else:
131 a1 = numpy.asarray(a)
132 if len(a1.shape) != 2:
133 raise ValueError("expected a 2-D array")
134 M, N = a1.shape
135 overwrite_a = overwrite_a or (_datacopied(a1, a))
137 if pivoting:
138 geqp3, = get_lapack_funcs(('geqp3',), (a1,))
139 qr, jpvt, tau = safecall(geqp3, "geqp3", a1, overwrite_a=overwrite_a)
140 jpvt -= 1 # geqp3 returns a 1-based index array, so subtract 1
141 else:
142 geqrf, = get_lapack_funcs(('geqrf',), (a1,))
143 qr, tau = safecall(geqrf, "geqrf", a1, lwork=lwork,
144 overwrite_a=overwrite_a)
146 if mode not in ['economic', 'raw'] or M < N:
147 R = numpy.triu(qr)
148 else:
149 R = numpy.triu(qr[:N, :])
151 if pivoting:
152 Rj = R, jpvt
153 else:
154 Rj = R,
156 if mode == 'r':
157 return Rj
158 elif mode == 'raw':
159 return ((qr, tau),) + Rj
161 gor_un_gqr, = get_lapack_funcs(('orgqr',), (qr,))
163 if M < N:
164 Q, = safecall(gor_un_gqr, "gorgqr/gungqr", qr[:, :M], tau,
165 lwork=lwork, overwrite_a=1)
166 elif mode == 'economic':
167 Q, = safecall(gor_un_gqr, "gorgqr/gungqr", qr, tau, lwork=lwork,
168 overwrite_a=1)
169 else:
170 t = qr.dtype.char
171 qqr = numpy.empty((M, M), dtype=t)
172 qqr[:, :N] = qr
173 Q, = safecall(gor_un_gqr, "gorgqr/gungqr", qqr, tau, lwork=lwork,
174 overwrite_a=1)
176 return (Q,) + Rj
179def qr_multiply(a, c, mode='right', pivoting=False, conjugate=False,
180 overwrite_a=False, overwrite_c=False):
181 """
182 Calculate the QR decomposition and multiply Q with a matrix.
184 Calculate the decomposition ``A = Q R`` where Q is unitary/orthogonal
185 and R upper triangular. Multiply Q with a vector or a matrix c.
187 Parameters
188 ----------
189 a : (M, N), array_like
190 Input array
191 c : array_like
192 Input array to be multiplied by ``q``.
193 mode : {'left', 'right'}, optional
194 ``Q @ c`` is returned if mode is 'left', ``c @ Q`` is returned if
195 mode is 'right'.
196 The shape of c must be appropriate for the matrix multiplications,
197 if mode is 'left', ``min(a.shape) == c.shape[0]``,
198 if mode is 'right', ``a.shape[0] == c.shape[1]``.
199 pivoting : bool, optional
200 Whether or not factorization should include pivoting for rank-revealing
201 qr decomposition, see the documentation of qr.
202 conjugate : bool, optional
203 Whether Q should be complex-conjugated. This might be faster
204 than explicit conjugation.
205 overwrite_a : bool, optional
206 Whether data in a is overwritten (may improve performance)
207 overwrite_c : bool, optional
208 Whether data in c is overwritten (may improve performance).
209 If this is used, c must be big enough to keep the result,
210 i.e. ``c.shape[0]`` = ``a.shape[0]`` if mode is 'left'.
212 Returns
213 -------
214 CQ : ndarray
215 The product of ``Q`` and ``c``.
216 R : (K, N), ndarray
217 R array of the resulting QR factorization where ``K = min(M, N)``.
218 P : (N,) ndarray
219 Integer pivot array. Only returned when ``pivoting=True``.
221 Raises
222 ------
223 LinAlgError
224 Raised if QR decomposition fails.
226 Notes
227 -----
228 This is an interface to the LAPACK routines ``?GEQRF``, ``?ORMQR``,
229 ``?UNMQR``, and ``?GEQP3``.
231 .. versionadded:: 0.11.0
233 Examples
234 --------
235 >>> import numpy as np
236 >>> from scipy.linalg import qr_multiply, qr
237 >>> A = np.array([[1, 3, 3], [2, 3, 2], [2, 3, 3], [1, 3, 2]])
238 >>> qc, r1, piv1 = qr_multiply(A, 2*np.eye(4), pivoting=1)
239 >>> qc
240 array([[-1., 1., -1.],
241 [-1., -1., 1.],
242 [-1., -1., -1.],
243 [-1., 1., 1.]])
244 >>> r1
245 array([[-6., -3., -5. ],
246 [ 0., -1., -1.11022302e-16],
247 [ 0., 0., -1. ]])
248 >>> piv1
249 array([1, 0, 2], dtype=int32)
250 >>> q2, r2, piv2 = qr(A, mode='economic', pivoting=1)
251 >>> np.allclose(2*q2 - qc, np.zeros((4, 3)))
252 True
254 """
255 if mode not in ['left', 'right']:
256 raise ValueError("Mode argument can only be 'left' or 'right' but "
257 "not '{}'".format(mode))
258 c = numpy.asarray_chkfinite(c)
259 if c.ndim < 2:
260 onedim = True
261 c = numpy.atleast_2d(c)
262 if mode == "left":
263 c = c.T
264 else:
265 onedim = False
267 a = numpy.atleast_2d(numpy.asarray(a)) # chkfinite done in qr
268 M, N = a.shape
270 if mode == 'left':
271 if c.shape[0] != min(M, N + overwrite_c*(M-N)):
272 raise ValueError('Array shapes are not compatible for Q @ c'
273 ' operation: {} vs {}'.format(a.shape, c.shape))
274 else:
275 if M != c.shape[1]:
276 raise ValueError('Array shapes are not compatible for c @ Q'
277 ' operation: {} vs {}'.format(c.shape, a.shape))
279 raw = qr(a, overwrite_a, None, "raw", pivoting)
280 Q, tau = raw[0]
282 gor_un_mqr, = get_lapack_funcs(('ormqr',), (Q,))
283 if gor_un_mqr.typecode in ('s', 'd'):
284 trans = "T"
285 else:
286 trans = "C"
288 Q = Q[:, :min(M, N)]
289 if M > N and mode == "left" and not overwrite_c:
290 if conjugate:
291 cc = numpy.zeros((c.shape[1], M), dtype=c.dtype, order="F")
292 cc[:, :N] = c.T
293 else:
294 cc = numpy.zeros((M, c.shape[1]), dtype=c.dtype, order="F")
295 cc[:N, :] = c
296 trans = "N"
297 if conjugate:
298 lr = "R"
299 else:
300 lr = "L"
301 overwrite_c = True
302 elif c.flags["C_CONTIGUOUS"] and trans == "T" or conjugate:
303 cc = c.T
304 if mode == "left":
305 lr = "R"
306 else:
307 lr = "L"
308 else:
309 trans = "N"
310 cc = c
311 if mode == "left":
312 lr = "L"
313 else:
314 lr = "R"
315 cQ, = safecall(gor_un_mqr, "gormqr/gunmqr", lr, trans, Q, tau, cc,
316 overwrite_c=overwrite_c)
317 if trans != "N":
318 cQ = cQ.T
319 if mode == "right":
320 cQ = cQ[:, :min(M, N)]
321 if onedim:
322 cQ = cQ.ravel()
324 return (cQ,) + raw[1:]
327def rq(a, overwrite_a=False, lwork=None, mode='full', check_finite=True):
328 """
329 Compute RQ decomposition of a matrix.
331 Calculate the decomposition ``A = R Q`` where Q is unitary/orthogonal
332 and R upper triangular.
334 Parameters
335 ----------
336 a : (M, N) array_like
337 Matrix to be decomposed
338 overwrite_a : bool, optional
339 Whether data in a is overwritten (may improve performance)
340 lwork : int, optional
341 Work array size, lwork >= a.shape[1]. If None or -1, an optimal size
342 is computed.
343 mode : {'full', 'r', 'economic'}, optional
344 Determines what information is to be returned: either both Q and R
345 ('full', default), only R ('r') or both Q and R but computed in
346 economy-size ('economic', see Notes).
347 check_finite : bool, optional
348 Whether to check that the input matrix contains only finite numbers.
349 Disabling may give a performance gain, but may result in problems
350 (crashes, non-termination) if the inputs do contain infinities or NaNs.
352 Returns
353 -------
354 R : float or complex ndarray
355 Of shape (M, N) or (M, K) for ``mode='economic'``. ``K = min(M, N)``.
356 Q : float or complex ndarray
357 Of shape (N, N) or (K, N) for ``mode='economic'``. Not returned
358 if ``mode='r'``.
360 Raises
361 ------
362 LinAlgError
363 If decomposition fails.
365 Notes
366 -----
367 This is an interface to the LAPACK routines sgerqf, dgerqf, cgerqf, zgerqf,
368 sorgrq, dorgrq, cungrq and zungrq.
370 If ``mode=economic``, the shapes of Q and R are (K, N) and (M, K) instead
371 of (N,N) and (M,N), with ``K=min(M,N)``.
373 Examples
374 --------
375 >>> import numpy as np
376 >>> from scipy import linalg
377 >>> rng = np.random.default_rng()
378 >>> a = rng.standard_normal((6, 9))
379 >>> r, q = linalg.rq(a)
380 >>> np.allclose(a, r @ q)
381 True
382 >>> r.shape, q.shape
383 ((6, 9), (9, 9))
384 >>> r2 = linalg.rq(a, mode='r')
385 >>> np.allclose(r, r2)
386 True
387 >>> r3, q3 = linalg.rq(a, mode='economic')
388 >>> r3.shape, q3.shape
389 ((6, 6), (6, 9))
391 """
392 if mode not in ['full', 'r', 'economic']:
393 raise ValueError(
394 "Mode argument should be one of ['full', 'r', 'economic']")
396 if check_finite:
397 a1 = numpy.asarray_chkfinite(a)
398 else:
399 a1 = numpy.asarray(a)
400 if len(a1.shape) != 2:
401 raise ValueError('expected matrix')
402 M, N = a1.shape
403 overwrite_a = overwrite_a or (_datacopied(a1, a))
405 gerqf, = get_lapack_funcs(('gerqf',), (a1,))
406 rq, tau = safecall(gerqf, 'gerqf', a1, lwork=lwork,
407 overwrite_a=overwrite_a)
408 if not mode == 'economic' or N < M:
409 R = numpy.triu(rq, N-M)
410 else:
411 R = numpy.triu(rq[-M:, -M:])
413 if mode == 'r':
414 return R
416 gor_un_grq, = get_lapack_funcs(('orgrq',), (rq,))
418 if N < M:
419 Q, = safecall(gor_un_grq, "gorgrq/gungrq", rq[-N:], tau, lwork=lwork,
420 overwrite_a=1)
421 elif mode == 'economic':
422 Q, = safecall(gor_un_grq, "gorgrq/gungrq", rq, tau, lwork=lwork,
423 overwrite_a=1)
424 else:
425 rq1 = numpy.empty((N, N), dtype=rq.dtype)
426 rq1[-M:] = rq
427 Q, = safecall(gor_un_grq, "gorgrq/gungrq", rq1, tau, lwork=lwork,
428 overwrite_a=1)
430 return R, Q