Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/linalg/decomp_qr.py : 6%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
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'``.
66 R : float or complex ndarray
67 Of shape (M, N), or (K, N) for ``mode='economic'``. ``K = min(M, N)``.
68 P : int ndarray
69 Of shape (N,) for ``pivoting=True``. Not returned if
70 ``pivoting=False``.
72 Raises
73 ------
74 LinAlgError
75 Raised if decomposition fails
77 Notes
78 -----
79 This is an interface to the LAPACK routines dgeqrf, zgeqrf,
80 dorgqr, zungqr, dgeqp3, and zgeqp3.
82 If ``mode=economic``, the shapes of Q and R are (M, K) and (K, N) instead
83 of (M,M) and (M,N), with ``K=min(M,N)``.
85 Examples
86 --------
87 >>> from scipy import linalg
88 >>> a = np.random.randn(9, 6)
90 >>> q, r = linalg.qr(a)
91 >>> np.allclose(a, np.dot(q, r))
92 True
93 >>> q.shape, r.shape
94 ((9, 9), (9, 6))
96 >>> r2 = linalg.qr(a, mode='r')
97 >>> np.allclose(r, r2)
98 True
100 >>> q3, r3 = linalg.qr(a, mode='economic')
101 >>> q3.shape, r3.shape
102 ((9, 6), (6, 6))
104 >>> q4, r4, p4 = linalg.qr(a, pivoting=True)
105 >>> d = np.abs(np.diag(r4))
106 >>> np.all(d[1:] <= d[:-1])
107 True
108 >>> np.allclose(a[:, p4], np.dot(q4, r4))
109 True
110 >>> q4.shape, r4.shape, p4.shape
111 ((9, 9), (9, 6), (6,))
113 >>> q5, r5, p5 = linalg.qr(a, mode='economic', pivoting=True)
114 >>> q5.shape, r5.shape, p5.shape
115 ((9, 6), (6, 6), (6,))
117 """
118 # 'qr' was the old default, equivalent to 'full'. Neither 'full' nor
119 # 'qr' are used below.
120 # 'raw' is used internally by qr_multiply
121 if mode not in ['full', 'qr', 'r', 'economic', 'raw']:
122 raise ValueError("Mode argument should be one of ['full', 'r',"
123 "'economic', 'raw']")
125 if check_finite:
126 a1 = numpy.asarray_chkfinite(a)
127 else:
128 a1 = numpy.asarray(a)
129 if len(a1.shape) != 2:
130 raise ValueError("expected a 2-D array")
131 M, N = a1.shape
132 overwrite_a = overwrite_a or (_datacopied(a1, a))
134 if pivoting:
135 geqp3, = get_lapack_funcs(('geqp3',), (a1,))
136 qr, jpvt, tau = safecall(geqp3, "geqp3", a1, overwrite_a=overwrite_a)
137 jpvt -= 1 # geqp3 returns a 1-based index array, so subtract 1
138 else:
139 geqrf, = get_lapack_funcs(('geqrf',), (a1,))
140 qr, tau = safecall(geqrf, "geqrf", a1, lwork=lwork,
141 overwrite_a=overwrite_a)
143 if mode not in ['economic', 'raw'] or M < N:
144 R = numpy.triu(qr)
145 else:
146 R = numpy.triu(qr[:N, :])
148 if pivoting:
149 Rj = R, jpvt
150 else:
151 Rj = R,
153 if mode == 'r':
154 return Rj
155 elif mode == 'raw':
156 return ((qr, tau),) + Rj
158 gor_un_gqr, = get_lapack_funcs(('orgqr',), (qr,))
160 if M < N:
161 Q, = safecall(gor_un_gqr, "gorgqr/gungqr", qr[:, :M], tau,
162 lwork=lwork, overwrite_a=1)
163 elif mode == 'economic':
164 Q, = safecall(gor_un_gqr, "gorgqr/gungqr", qr, tau, lwork=lwork,
165 overwrite_a=1)
166 else:
167 t = qr.dtype.char
168 qqr = numpy.empty((M, M), dtype=t)
169 qqr[:, :N] = qr
170 Q, = safecall(gor_un_gqr, "gorgqr/gungqr", qqr, tau, lwork=lwork,
171 overwrite_a=1)
173 return (Q,) + Rj
176def qr_multiply(a, c, mode='right', pivoting=False, conjugate=False,
177 overwrite_a=False, overwrite_c=False):
178 """
179 Calculate the QR decomposition and multiply Q with a matrix.
181 Calculate the decomposition ``A = Q R`` where Q is unitary/orthogonal
182 and R upper triangular. Multiply Q with a vector or a matrix c.
184 Parameters
185 ----------
186 a : (M, N), array_like
187 Input array
188 c : array_like
189 Input array to be multiplied by ``q``.
190 mode : {'left', 'right'}, optional
191 ``Q @ c`` is returned if mode is 'left', ``c @ Q`` is returned if
192 mode is 'right'.
193 The shape of c must be appropriate for the matrix multiplications,
194 if mode is 'left', ``min(a.shape) == c.shape[0]``,
195 if mode is 'right', ``a.shape[0] == c.shape[1]``.
196 pivoting : bool, optional
197 Whether or not factorization should include pivoting for rank-revealing
198 qr decomposition, see the documentation of qr.
199 conjugate : bool, optional
200 Whether Q should be complex-conjugated. This might be faster
201 than explicit conjugation.
202 overwrite_a : bool, optional
203 Whether data in a is overwritten (may improve performance)
204 overwrite_c : bool, optional
205 Whether data in c is overwritten (may improve performance).
206 If this is used, c must be big enough to keep the result,
207 i.e. ``c.shape[0]`` = ``a.shape[0]`` if mode is 'left'.
209 Returns
210 -------
211 CQ : ndarray
212 The product of ``Q`` and ``c``.
213 R : (K, N), ndarray
214 R array of the resulting QR factorization where ``K = min(M, N)``.
215 P : (N,) ndarray
216 Integer pivot array. Only returned when ``pivoting=True``.
218 Raises
219 ------
220 LinAlgError
221 Raised if QR decomposition fails.
223 Notes
224 -----
225 This is an interface to the LAPACK routines ``?GEQRF``, ``?ORMQR``,
226 ``?UNMQR``, and ``?GEQP3``.
228 .. versionadded:: 0.11.0
230 Examples
231 --------
232 >>> from scipy.linalg import qr_multiply, qr
233 >>> A = np.array([[1, 3, 3], [2, 3, 2], [2, 3, 3], [1, 3, 2]])
234 >>> qc, r1, piv1 = qr_multiply(A, 2*np.eye(4), pivoting=1)
235 >>> qc
236 array([[-1., 1., -1.],
237 [-1., -1., 1.],
238 [-1., -1., -1.],
239 [-1., 1., 1.]])
240 >>> r1
241 array([[-6., -3., -5. ],
242 [ 0., -1., -1.11022302e-16],
243 [ 0., 0., -1. ]])
244 >>> piv1
245 array([1, 0, 2], dtype=int32)
246 >>> q2, r2, piv2 = qr(A, mode='economic', pivoting=1)
247 >>> np.allclose(2*q2 - qc, np.zeros((4, 3)))
248 True
250 """
251 if mode not in ['left', 'right']:
252 raise ValueError("Mode argument can only be 'left' or 'right' but "
253 "not '{}'".format(mode))
254 c = numpy.asarray_chkfinite(c)
255 if c.ndim < 2:
256 onedim = True
257 c = numpy.atleast_2d(c)
258 if mode == "left":
259 c = c.T
260 else:
261 onedim = False
263 a = numpy.atleast_2d(numpy.asarray(a)) # chkfinite done in qr
264 M, N = a.shape
266 if mode == 'left':
267 if c.shape[0] != min(M, N + overwrite_c*(M-N)):
268 raise ValueError('Array shapes are not compatible for Q @ c'
269 ' operation: {} vs {}'.format(a.shape, c.shape))
270 else:
271 if M != c.shape[1]:
272 raise ValueError('Array shapes are not compatible for c @ Q'
273 ' operation: {} vs {}'.format(c.shape, a.shape))
275 raw = qr(a, overwrite_a, None, "raw", pivoting)
276 Q, tau = raw[0]
278 gor_un_mqr, = get_lapack_funcs(('ormqr',), (Q,))
279 if gor_un_mqr.typecode in ('s', 'd'):
280 trans = "T"
281 else:
282 trans = "C"
284 Q = Q[:, :min(M, N)]
285 if M > N and mode == "left" and not overwrite_c:
286 if conjugate:
287 cc = numpy.zeros((c.shape[1], M), dtype=c.dtype, order="F")
288 cc[:, :N] = c.T
289 else:
290 cc = numpy.zeros((M, c.shape[1]), dtype=c.dtype, order="F")
291 cc[:N, :] = c
292 trans = "N"
293 if conjugate:
294 lr = "R"
295 else:
296 lr = "L"
297 overwrite_c = True
298 elif c.flags["C_CONTIGUOUS"] and trans == "T" or conjugate:
299 cc = c.T
300 if mode == "left":
301 lr = "R"
302 else:
303 lr = "L"
304 else:
305 trans = "N"
306 cc = c
307 if mode == "left":
308 lr = "L"
309 else:
310 lr = "R"
311 cQ, = safecall(gor_un_mqr, "gormqr/gunmqr", lr, trans, Q, tau, cc,
312 overwrite_c=overwrite_c)
313 if trans != "N":
314 cQ = cQ.T
315 if mode == "right":
316 cQ = cQ[:, :min(M, N)]
317 if onedim:
318 cQ = cQ.ravel()
320 return (cQ,) + raw[1:]
323def rq(a, overwrite_a=False, lwork=None, mode='full', check_finite=True):
324 """
325 Compute RQ decomposition of a matrix.
327 Calculate the decomposition ``A = R Q`` where Q is unitary/orthogonal
328 and R upper triangular.
330 Parameters
331 ----------
332 a : (M, N) array_like
333 Matrix to be decomposed
334 overwrite_a : bool, optional
335 Whether data in a is overwritten (may improve performance)
336 lwork : int, optional
337 Work array size, lwork >= a.shape[1]. If None or -1, an optimal size
338 is computed.
339 mode : {'full', 'r', 'economic'}, optional
340 Determines what information is to be returned: either both Q and R
341 ('full', default), only R ('r') or both Q and R but computed in
342 economy-size ('economic', see Notes).
343 check_finite : bool, optional
344 Whether to check that the input matrix contains only finite numbers.
345 Disabling may give a performance gain, but may result in problems
346 (crashes, non-termination) if the inputs do contain infinities or NaNs.
348 Returns
349 -------
350 R : float or complex ndarray
351 Of shape (M, N) or (M, K) for ``mode='economic'``. ``K = min(M, N)``.
352 Q : float or complex ndarray
353 Of shape (N, N) or (K, N) for ``mode='economic'``. Not returned
354 if ``mode='r'``.
356 Raises
357 ------
358 LinAlgError
359 If decomposition fails.
361 Notes
362 -----
363 This is an interface to the LAPACK routines sgerqf, dgerqf, cgerqf, zgerqf,
364 sorgrq, dorgrq, cungrq and zungrq.
366 If ``mode=economic``, the shapes of Q and R are (K, N) and (M, K) instead
367 of (N,N) and (M,N), with ``K=min(M,N)``.
369 Examples
370 --------
371 >>> from scipy import linalg
372 >>> a = np.random.randn(6, 9)
373 >>> r, q = linalg.rq(a)
374 >>> np.allclose(a, r @ q)
375 True
376 >>> r.shape, q.shape
377 ((6, 9), (9, 9))
378 >>> r2 = linalg.rq(a, mode='r')
379 >>> np.allclose(r, r2)
380 True
381 >>> r3, q3 = linalg.rq(a, mode='economic')
382 >>> r3.shape, q3.shape
383 ((6, 6), (6, 9))
385 """
386 if mode not in ['full', 'r', 'economic']:
387 raise ValueError(
388 "Mode argument should be one of ['full', 'r', 'economic']")
390 if check_finite:
391 a1 = numpy.asarray_chkfinite(a)
392 else:
393 a1 = numpy.asarray(a)
394 if len(a1.shape) != 2:
395 raise ValueError('expected matrix')
396 M, N = a1.shape
397 overwrite_a = overwrite_a or (_datacopied(a1, a))
399 gerqf, = get_lapack_funcs(('gerqf',), (a1,))
400 rq, tau = safecall(gerqf, 'gerqf', a1, lwork=lwork,
401 overwrite_a=overwrite_a)
402 if not mode == 'economic' or N < M:
403 R = numpy.triu(rq, N-M)
404 else:
405 R = numpy.triu(rq[-M:, -M:])
407 if mode == 'r':
408 return R
410 gor_un_grq, = get_lapack_funcs(('orgrq',), (rq,))
412 if N < M:
413 Q, = safecall(gor_un_grq, "gorgrq/gungrq", rq[-N:], tau, lwork=lwork,
414 overwrite_a=1)
415 elif mode == 'economic':
416 Q, = safecall(gor_un_grq, "gorgrq/gungrq", rq, tau, lwork=lwork,
417 overwrite_a=1)
418 else:
419 rq1 = numpy.empty((N, N), dtype=rq.dtype)
420 rq1[-M:] = rq
421 Q, = safecall(gor_un_grq, "gorgrq/gungrq", rq1, tau, lwork=lwork,
422 overwrite_a=1)
424 return R, Q