Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/linalg/_decomp

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 >>> import numpy as np

88 >>> from scipy import linalg

89 >>> rng = np.random.default_rng()

90 >>> a = rng.standard_normal((9, 6))

92 >>> q, r = linalg.qr(a)

93 >>> np.allclose(a, np.dot(q, r))

94 True

95 >>> q.shape, r.shape

96 ((9, 9), (9, 6))

98 >>> r2 = linalg.qr(a, mode='r')

99 >>> np.allclose(r, r2)

100 True

101

102 >>> q3, r3 = linalg.qr(a, mode='economic')

103 >>> q3.shape, r3.shape

104 ((9, 6), (6, 6))

105

106 >>> q4, r4, p4 = linalg.qr(a, pivoting=True)

107 >>> d = np.abs(np.diag(r4))

108 >>> np.all(d[1:] <= d[:-1])

109 True

110 >>> np.allclose(a[:, p4], np.dot(q4, r4))

111 True

112 >>> q4.shape, r4.shape, p4.shape

113 ((9, 9), (9, 6), (6,))

114

115 >>> q5, r5, p5 = linalg.qr(a, mode='economic', pivoting=True)

116 >>> q5.shape, r5.shape, p5.shape

117 ((9, 6), (6, 6), (6,))

118

119 """

120 # 'qr' was the old default, equivalent to 'full'. Neither 'full' nor

121 # 'qr' are used below.

122 # 'raw' is used internally by qr_multiply

123 if mode not in ['full', 'qr', 'r', 'economic', 'raw']:

124 raise ValueError("Mode argument should be one of ['full', 'r',"

125 "'economic', 'raw']")

126

127 if check_finite:

128 a1 = numpy.asarray_chkfinite(a)

129 else:

130 a1 = numpy.asarray(a)

131 if len(a1.shape) != 2:

132 raise ValueError("expected a 2-D array")

133 M, N = a1.shape

134 overwrite_a = overwrite_a or (_datacopied(a1, a))

135

136 if pivoting:

137 geqp3, = get_lapack_funcs(('geqp3',), (a1,))

138 qr, jpvt, tau = safecall(geqp3, "geqp3", a1, overwrite_a=overwrite_a)

139 jpvt -= 1 # geqp3 returns a 1-based index array, so subtract 1

140 else:

141 geqrf, = get_lapack_funcs(('geqrf',), (a1,))

142 qr, tau = safecall(geqrf, "geqrf", a1, lwork=lwork,

143 overwrite_a=overwrite_a)

144

145 if mode not in ['economic', 'raw'] or M < N:

146 R = numpy.triu(qr)

147 else:

148 R = numpy.triu(qr[:N, :])

149

150 if pivoting:

151 Rj = R, jpvt

152 else:

153 Rj = R,

154

155 if mode == 'r':

156 return Rj

157 elif mode == 'raw':

158 return ((qr, tau),) + Rj

159

160 gor_un_gqr, = get_lapack_funcs(('orgqr',), (qr,))

161

162 if M < N:

163 Q, = safecall(gor_un_gqr, "gorgqr/gungqr", qr[:, :M], tau,

164 lwork=lwork, overwrite_a=1)

165 elif mode == 'economic':

166 Q, = safecall(gor_un_gqr, "gorgqr/gungqr", qr, tau, lwork=lwork,

167 overwrite_a=1)

168 else:

169 t = qr.dtype.char

170 qqr = numpy.empty((M, M), dtype=t)

171 qqr[:, :N] = qr

172 Q, = safecall(gor_un_gqr, "gorgqr/gungqr", qqr, tau, lwork=lwork,

173 overwrite_a=1)

174

175 return (Q,) + Rj

176

177

178def qr_multiply(a, c, mode='right', pivoting=False, conjugate=False,

179 overwrite_a=False, overwrite_c=False):

180 """

181 Calculate the QR decomposition and multiply Q with a matrix.

182

183 Calculate the decomposition ``A = Q R`` where Q is unitary/orthogonal

184 and R upper triangular. Multiply Q with a vector or a matrix c.

185

186 Parameters

187 ----------

188 a : (M, N), array_like

189 Input array

190 c : array_like

191 Input array to be multiplied by ``q``.

192 mode : {'left', 'right'}, optional

193 ``Q @ c`` is returned if mode is 'left', ``c @ Q`` is returned if

194 mode is 'right'.

195 The shape of c must be appropriate for the matrix multiplications,

196 if mode is 'left', ``min(a.shape) == c.shape[0]``,

197 if mode is 'right', ``a.shape[0] == c.shape[1]``.

198 pivoting : bool, optional

199 Whether or not factorization should include pivoting for rank-revealing

200 qr decomposition, see the documentation of qr.

201 conjugate : bool, optional

202 Whether Q should be complex-conjugated. This might be faster

203 than explicit conjugation.

204 overwrite_a : bool, optional

205 Whether data in a is overwritten (may improve performance)

206 overwrite_c : bool, optional

207 Whether data in c is overwritten (may improve performance).

208 If this is used, c must be big enough to keep the result,

209 i.e. ``c.shape[0]`` = ``a.shape[0]`` if mode is 'left'.

210

211 Returns

212 -------

213 CQ : ndarray

214 The product of ``Q`` and ``c``.

215 R : (K, N), ndarray

216 R array of the resulting QR factorization where ``K = min(M, N)``.

217 P : (N,) ndarray

218 Integer pivot array. Only returned when ``pivoting=True``.

219

220 Raises

221 ------

222 LinAlgError

223 Raised if QR decomposition fails.

224

225 Notes

226 -----

227 This is an interface to the LAPACK routines ``?GEQRF``, ``?ORMQR``,

228 ``?UNMQR``, and ``?GEQP3``.

229

230 .. versionadded:: 0.11.0

231

232 Examples

233 --------

234 >>> import numpy as np

235 >>> from scipy.linalg import qr_multiply, qr

236 >>> A = np.array([[1, 3, 3], [2, 3, 2], [2, 3, 3], [1, 3, 2]])

237 >>> qc, r1, piv1 = qr_multiply(A, 2*np.eye(4), pivoting=1)

238 >>> qc

239 array([[-1., 1., -1.],

240 [-1., -1., 1.],

241 [-1., -1., -1.],

242 [-1., 1., 1.]])

243 >>> r1

244 array([[-6., -3., -5. ],

245 [ 0., -1., -1.11022302e-16],

246 [ 0., 0., -1. ]])

247 >>> piv1

248 array([1, 0, 2], dtype=int32)

249 >>> q2, r2, piv2 = qr(A, mode='economic', pivoting=1)

250 >>> np.allclose(2*q2 - qc, np.zeros((4, 3)))

251 True

252

253 """

254 if mode not in ['left', 'right']:

255 raise ValueError("Mode argument can only be 'left' or 'right' but "

256 "not '{}'".format(mode))

257 c = numpy.asarray_chkfinite(c)

258 if c.ndim < 2:

259 onedim = True

260 c = numpy.atleast_2d(c)

261 if mode == "left":

262 c = c.T

263 else:

264 onedim = False

265

266 a = numpy.atleast_2d(numpy.asarray(a)) # chkfinite done in qr

267 M, N = a.shape

268

269 if mode == 'left':

270 if c.shape[0] != min(M, N + overwrite_c*(M-N)):

271 raise ValueError('Array shapes are not compatible for Q @ c'

272 ' operation: {} vs {}'.format(a.shape, c.shape))

273 else:

274 if M != c.shape[1]:

275 raise ValueError('Array shapes are not compatible for c @ Q'

276 ' operation: {} vs {}'.format(c.shape, a.shape))

277

278 raw = qr(a, overwrite_a, None, "raw", pivoting)

279 Q, tau = raw[0]

280

281 gor_un_mqr, = get_lapack_funcs(('ormqr',), (Q,))

282 if gor_un_mqr.typecode in ('s', 'd'):

283 trans = "T"

284 else:

285 trans = "C"

286

287 Q = Q[:, :min(M, N)]

288 if M > N and mode == "left" and not overwrite_c:

289 if conjugate:

290 cc = numpy.zeros((c.shape[1], M), dtype=c.dtype, order="F")

291 cc[:, :N] = c.T

292 else:

293 cc = numpy.zeros((M, c.shape[1]), dtype=c.dtype, order="F")

294 cc[:N, :] = c

295 trans = "N"

296 if conjugate:

297 lr = "R"

298 else:

299 lr = "L"

300 overwrite_c = True

301 elif c.flags["C_CONTIGUOUS"] and trans == "T" or conjugate:

302 cc = c.T

303 if mode == "left":

304 lr = "R"

305 else:

306 lr = "L"

307 else:

308 trans = "N"

309 cc = c

310 if mode == "left":

311 lr = "L"

312 else:

313 lr = "R"

314 cQ, = safecall(gor_un_mqr, "gormqr/gunmqr", lr, trans, Q, tau, cc,

315 overwrite_c=overwrite_c)

316 if trans != "N":

317 cQ = cQ.T

318 if mode == "right":

319 cQ = cQ[:, :min(M, N)]

320 if onedim:

321 cQ = cQ.ravel()

322

323 return (cQ,) + raw[1:]

324

325

326def rq(a, overwrite_a=False, lwork=None, mode='full', check_finite=True):

327 """

328 Compute RQ decomposition of a matrix.

329

330 Calculate the decomposition ``A = R Q`` where Q is unitary/orthogonal

331 and R upper triangular.

332

333 Parameters

334 ----------

335 a : (M, N) array_like

336 Matrix to be decomposed

337 overwrite_a : bool, optional

338 Whether data in a is overwritten (may improve performance)

339 lwork : int, optional

340 Work array size, lwork >= a.shape[1]. If None or -1, an optimal size

341 is computed.

342 mode : {'full', 'r', 'economic'}, optional

343 Determines what information is to be returned: either both Q and R

344 ('full', default), only R ('r') or both Q and R but computed in

345 economy-size ('economic', see Notes).

346 check_finite : bool, optional

347 Whether to check that the input matrix contains only finite numbers.

348 Disabling may give a performance gain, but may result in problems

349 (crashes, non-termination) if the inputs do contain infinities or NaNs.

350

351 Returns

352 -------

353 R : float or complex ndarray

354 Of shape (M, N) or (M, K) for ``mode='economic'``. ``K = min(M, N)``.

355 Q : float or complex ndarray

356 Of shape (N, N) or (K, N) for ``mode='economic'``. Not returned

357 if ``mode='r'``.

358

359 Raises

360 ------

361 LinAlgError

362 If decomposition fails.

363

364 Notes

365 -----

366 This is an interface to the LAPACK routines sgerqf, dgerqf, cgerqf, zgerqf,

367 sorgrq, dorgrq, cungrq and zungrq.

368

369 If ``mode=economic``, the shapes of Q and R are (K, N) and (M, K) instead

370 of (N,N) and (M,N), with ``K=min(M,N)``.

371

372 Examples

373 --------

374 >>> import numpy as np

375 >>> from scipy import linalg

376 >>> rng = np.random.default_rng()

377 >>> a = rng.standard_normal((6, 9))

378 >>> r, q = linalg.rq(a)

379 >>> np.allclose(a, r @ q)

380 True

381 >>> r.shape, q.shape

382 ((6, 9), (9, 9))

383 >>> r2 = linalg.rq(a, mode='r')

384 >>> np.allclose(r, r2)

385 True

386 >>> r3, q3 = linalg.rq(a, mode='economic')

387 >>> r3.shape, q3.shape

388 ((6, 6), (6, 9))

389

390 """

391 if mode not in ['full', 'r', 'economic']:

392 raise ValueError(

393 "Mode argument should be one of ['full', 'r', 'economic']")

394

395 if check_finite:

396 a1 = numpy.asarray_chkfinite(a)

397 else:

398 a1 = numpy.asarray(a)

399 if len(a1.shape) != 2:

400 raise ValueError('expected matrix')

401 M, N = a1.shape

402 overwrite_a = overwrite_a or (_datacopied(a1, a))

403

404 gerqf, = get_lapack_funcs(('gerqf',), (a1,))

405 rq, tau = safecall(gerqf, 'gerqf', a1, lwork=lwork,

406 overwrite_a=overwrite_a)

407 if not mode == 'economic' or N < M:

408 R = numpy.triu(rq, N-M)

409 else:

410 R = numpy.triu(rq[-M:, -M:])

411

412 if mode == 'r':

413 return R

414

415 gor_un_grq, = get_lapack_funcs(('orgrq',), (rq,))

416

417 if N < M:

418 Q, = safecall(gor_un_grq, "gorgrq/gungrq", rq[-N:], tau, lwork=lwork,

419 overwrite_a=1)

420 elif mode == 'economic':

421 Q, = safecall(gor_un_grq, "gorgrq/gungrq", rq, tau, lwork=lwork,

422 overwrite_a=1)

423 else:

424 rq1 = numpy.empty((N, N), dtype=rq.dtype)

425 rq1[-M:] = rq

426 Q, = safecall(gor_un_grq, "gorgrq/gungrq", rq1, tau, lwork=lwork,

427 overwrite_a=1)

428

429 return R, Q

Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/linalg/_decomp_qr.py: 6%

129 statements