Coverage for /usr/lib/python3/dist-packages/mpmath/matrices/eigen_symmetric.py: 3%
824 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1#!/usr/bin/python
2# -*- coding: utf-8 -*-
4##################################################################################################
5# module for the symmetric eigenvalue problem
6# Copyright 2013 Timo Hartmann (thartmann15 at gmail.com)
7#
8# todo:
9# - implement balancing
10#
11##################################################################################################
13"""
14The symmetric eigenvalue problem.
15---------------------------------
17This file contains routines for the symmetric eigenvalue problem.
19high level routines:
21 eigsy : real symmetric (ordinary) eigenvalue problem
22 eighe : complex hermitian (ordinary) eigenvalue problem
23 eigh : unified interface for eigsy and eighe
24 svd_r : singular value decomposition for real matrices
25 svd_c : singular value decomposition for complex matrices
26 svd : unified interface for svd_r and svd_c
29low level routines:
31 r_sy_tridiag : reduction of real symmetric matrix to real symmetric tridiagonal matrix
32 c_he_tridiag_0 : reduction of complex hermitian matrix to real symmetric tridiagonal matrix
33 c_he_tridiag_1 : auxiliary routine to c_he_tridiag_0
34 c_he_tridiag_2 : auxiliary routine to c_he_tridiag_0
35 tridiag_eigen : solves the real symmetric tridiagonal matrix eigenvalue problem
36 svd_r_raw : raw singular value decomposition for real matrices
37 svd_c_raw : raw singular value decomposition for complex matrices
38"""
40from ..libmp.backend import xrange
41from .eigen import defun
44def r_sy_tridiag(ctx, A, D, E, calc_ev = True):
45 """
46 This routine transforms a real symmetric matrix A to a real symmetric
47 tridiagonal matrix T using an orthogonal similarity transformation:
48 Q' * A * Q = T (here ' denotes the matrix transpose).
49 The orthogonal matrix Q is build up from Householder reflectors.
51 parameters:
52 A (input/output) On input, A contains the real symmetric matrix of
53 dimension (n,n). On output, if calc_ev is true, A contains the
54 orthogonal matrix Q, otherwise A is destroyed.
56 D (output) real array of length n, contains the diagonal elements
57 of the tridiagonal matrix
59 E (output) real array of length n, contains the offdiagonal elements
60 of the tridiagonal matrix in E[0:(n-1)] where is the dimension of
61 the matrix A. E[n-1] is undefined.
63 calc_ev (input) If calc_ev is true, this routine explicitly calculates the
64 orthogonal matrix Q which is then returned in A. If calc_ev is
65 false, Q is not explicitly calculated resulting in a shorter run time.
67 This routine is a python translation of the fortran routine tred2.f in the
68 software library EISPACK (see netlib.org) which itself is based on the algol
69 procedure tred2 described in:
70 - Num. Math. 11, p.181-195 (1968) by Martin, Reinsch and Wilkonson
71 - Handbook for auto. comp., Vol II, Linear Algebra, p.212-226 (1971)
73 For a good introduction to Householder reflections, see also
74 Stoer, Bulirsch - Introduction to Numerical Analysis.
75 """
77 # note : the vector v of the i-th houshoulder reflector is stored in a[(i+1):,i]
78 # whereas v/<v,v> is stored in a[i,(i+1):]
80 n = A.rows
81 for i in xrange(n - 1, 0, -1):
82 # scale the vector
84 scale = 0
85 for k in xrange(0, i):
86 scale += abs(A[k,i])
88 scale_inv = 0
89 if scale != 0:
90 scale_inv = 1/scale
92 # sadly there are floating point numbers not equal to zero whose reciprocal is infinity
94 if i == 1 or scale == 0 or ctx.isinf(scale_inv):
95 E[i] = A[i-1,i] # nothing to do
96 D[i] = 0
97 continue
99 # calculate parameters for housholder transformation
101 H = 0
102 for k in xrange(0, i):
103 A[k,i] *= scale_inv
104 H += A[k,i] * A[k,i]
106 F = A[i-1,i]
107 G = ctx.sqrt(H)
108 if F > 0:
109 G = -G
110 E[i] = scale * G
111 H -= F * G
112 A[i-1,i] = F - G
113 F = 0
115 # apply housholder transformation
117 for j in xrange(0, i):
118 if calc_ev:
119 A[i,j] = A[j,i] / H
121 G = 0 # calculate A*U
122 for k in xrange(0, j + 1):
123 G += A[k,j] * A[k,i]
124 for k in xrange(j + 1, i):
125 G += A[j,k] * A[k,i]
127 E[j] = G / H # calculate P
128 F += E[j] * A[j,i]
130 HH = F / (2 * H)
132 for j in xrange(0, i): # calculate reduced A
133 F = A[j,i]
134 G = E[j] - HH * F # calculate Q
135 E[j] = G
137 for k in xrange(0, j + 1):
138 A[k,j] -= F * E[k] + G * A[k,i]
140 D[i] = H
142 for i in xrange(1, n): # better for compatibility
143 E[i-1] = E[i]
144 E[n-1] = 0
146 if calc_ev:
147 D[0] = 0
148 for i in xrange(0, n):
149 if D[i] != 0:
150 for j in xrange(0, i): # accumulate transformation matrices
151 G = 0
152 for k in xrange(0, i):
153 G += A[i,k] * A[k,j]
154 for k in xrange(0, i):
155 A[k,j] -= G * A[k,i]
157 D[i] = A[i,i]
158 A[i,i] = 1
160 for j in xrange(0, i):
161 A[j,i] = A[i,j] = 0
162 else:
163 for i in xrange(0, n):
164 D[i] = A[i,i]
170def c_he_tridiag_0(ctx, A, D, E, T):
171 """
172 This routine transforms a complex hermitian matrix A to a real symmetric
173 tridiagonal matrix T using an unitary similarity transformation:
174 Q' * A * Q = T (here ' denotes the hermitian matrix transpose,
175 i.e. transposition und conjugation).
176 The unitary matrix Q is build up from Householder reflectors and
177 an unitary diagonal matrix.
179 parameters:
180 A (input/output) On input, A contains the complex hermitian matrix
181 of dimension (n,n). On output, A contains the unitary matrix Q
182 in compressed form.
184 D (output) real array of length n, contains the diagonal elements
185 of the tridiagonal matrix.
187 E (output) real array of length n, contains the offdiagonal elements
188 of the tridiagonal matrix in E[0:(n-1)] where is the dimension of
189 the matrix A. E[n-1] is undefined.
191 T (output) complex array of length n, contains a unitary diagonal
192 matrix.
194 This routine is a python translation (in slightly modified form) of the fortran
195 routine htridi.f in the software library EISPACK (see netlib.org) which itself
196 is a complex version of the algol procedure tred1 described in:
197 - Num. Math. 11, p.181-195 (1968) by Martin, Reinsch and Wilkonson
198 - Handbook for auto. comp., Vol II, Linear Algebra, p.212-226 (1971)
200 For a good introduction to Householder reflections, see also
201 Stoer, Bulirsch - Introduction to Numerical Analysis.
202 """
204 n = A.rows
205 T[n-1] = 1
206 for i in xrange(n - 1, 0, -1):
208 # scale the vector
210 scale = 0
211 for k in xrange(0, i):
212 scale += abs(ctx.re(A[k,i])) + abs(ctx.im(A[k,i]))
214 scale_inv = 0
215 if scale != 0:
216 scale_inv = 1 / scale
218 # sadly there are floating point numbers not equal to zero whose reciprocal is infinity
220 if scale == 0 or ctx.isinf(scale_inv):
221 E[i] = 0
222 D[i] = 0
223 T[i-1] = 1
224 continue
226 if i == 1:
227 F = A[i-1,i]
228 f = abs(F)
229 E[i] = f
230 D[i] = 0
231 if f != 0:
232 T[i-1] = T[i] * F / f
233 else:
234 T[i-1] = T[i]
235 continue
237 # calculate parameters for housholder transformation
239 H = 0
240 for k in xrange(0, i):
241 A[k,i] *= scale_inv
242 rr = ctx.re(A[k,i])
243 ii = ctx.im(A[k,i])
244 H += rr * rr + ii * ii
246 F = A[i-1,i]
247 f = abs(F)
248 G = ctx.sqrt(H)
249 H += G * f
250 E[i] = scale * G
251 if f != 0:
252 F = F / f
253 TZ = - T[i] * F # T[i-1]=-T[i]*F, but we need T[i-1] as temporary storage
254 G *= F
255 else:
256 TZ = -T[i] # T[i-1]=-T[i]
257 A[i-1,i] += G
258 F = 0
260 # apply housholder transformation
262 for j in xrange(0, i):
263 A[i,j] = A[j,i] / H
265 G = 0 # calculate A*U
266 for k in xrange(0, j + 1):
267 G += ctx.conj(A[k,j]) * A[k,i]
268 for k in xrange(j + 1, i):
269 G += A[j,k] * A[k,i]
271 T[j] = G / H # calculate P
272 F += ctx.conj(T[j]) * A[j,i]
274 HH = F / (2 * H)
276 for j in xrange(0, i): # calculate reduced A
277 F = A[j,i]
278 G = T[j] - HH * F # calculate Q
279 T[j] = G
281 for k in xrange(0, j + 1):
282 A[k,j] -= ctx.conj(F) * T[k] + ctx.conj(G) * A[k,i]
283 # as we use the lower left part for storage
284 # we have to use the transpose of the normal formula
286 T[i-1] = TZ
287 D[i] = H
289 for i in xrange(1, n): # better for compatibility
290 E[i-1] = E[i]
291 E[n-1] = 0
293 D[0] = 0
294 for i in xrange(0, n):
295 zw = D[i]
296 D[i] = ctx.re(A[i,i])
297 A[i,i] = zw
305def c_he_tridiag_1(ctx, A, T):
306 """
307 This routine forms the unitary matrix Q described in c_he_tridiag_0.
309 parameters:
310 A (input/output) On input, A is the same matrix as delivered by
311 c_he_tridiag_0. On output, A is set to Q.
313 T (input) On input, T is the same array as delivered by c_he_tridiag_0.
315 """
317 n = A.rows
319 for i in xrange(0, n):
320 if A[i,i] != 0:
321 for j in xrange(0, i):
322 G = 0
323 for k in xrange(0, i):
324 G += ctx.conj(A[i,k]) * A[k,j]
325 for k in xrange(0, i):
326 A[k,j] -= G * A[k,i]
328 A[i,i] = 1
330 for j in xrange(0, i):
331 A[j,i] = A[i,j] = 0
333 for i in xrange(0, n):
334 for k in xrange(0, n):
335 A[i,k] *= T[k]
340def c_he_tridiag_2(ctx, A, T, B):
341 """
342 This routine applied the unitary matrix Q described in c_he_tridiag_0
343 onto the the matrix B, i.e. it forms Q*B.
345 parameters:
346 A (input) On input, A is the same matrix as delivered by c_he_tridiag_0.
348 T (input) On input, T is the same array as delivered by c_he_tridiag_0.
350 B (input/output) On input, B is a complex matrix. On output B is replaced
351 by Q*B.
353 This routine is a python translation of the fortran routine htribk.f in the
354 software library EISPACK (see netlib.org). See c_he_tridiag_0 for more
355 references.
356 """
358 n = A.rows
360 for i in xrange(0, n):
361 for k in xrange(0, n):
362 B[k,i] *= T[k]
364 for i in xrange(0, n):
365 if A[i,i] != 0:
366 for j in xrange(0, n):
367 G = 0
368 for k in xrange(0, i):
369 G += ctx.conj(A[i,k]) * B[k,j]
370 for k in xrange(0, i):
371 B[k,j] -= G * A[k,i]
377def tridiag_eigen(ctx, d, e, z = False):
378 """
379 This subroutine find the eigenvalues and the first components of the
380 eigenvectors of a real symmetric tridiagonal matrix using the implicit
381 QL method.
383 parameters:
385 d (input/output) real array of length n. on input, d contains the diagonal
386 elements of the input matrix. on output, d contains the eigenvalues in
387 ascending order.
389 e (input) real array of length n. on input, e contains the offdiagonal
390 elements of the input matrix in e[0:(n-1)]. On output, e has been
391 destroyed.
393 z (input/output) If z is equal to False, no eigenvectors will be computed.
394 Otherwise on input z should have the format z[0:m,0:n] (i.e. a real or
395 complex matrix of dimension (m,n) ). On output this matrix will be
396 multiplied by the matrix of the eigenvectors (i.e. the columns of this
397 matrix are the eigenvectors): z --> z*EV
398 That means if z[i,j]={1 if j==j; 0 otherwise} on input, then on output
399 z will contain the first m components of the eigenvectors. That means
400 if m is equal to n, the i-th eigenvector will be z[:,i].
402 This routine is a python translation (in slightly modified form) of the
403 fortran routine imtql2.f in the software library EISPACK (see netlib.org)
404 which itself is based on the algol procudure imtql2 desribed in:
405 - num. math. 12, p. 377-383(1968) by matrin and wilkinson
406 - modified in num. math. 15, p. 450(1970) by dubrulle
407 - handbook for auto. comp., vol. II-linear algebra, p. 241-248 (1971)
408 See also the routine gaussq.f in netlog.org or acm algorithm 726.
409 """
411 n = len(d)
412 e[n-1] = 0
413 iterlim = 2 * ctx.dps
415 for l in xrange(n):
416 j = 0
417 while 1:
418 m = l
419 while 1:
420 # look for a small subdiagonal element
421 if m + 1 == n:
422 break
423 if abs(e[m]) <= ctx.eps * (abs(d[m]) + abs(d[m + 1])):
424 break
425 m = m + 1
426 if m == l:
427 break
429 if j >= iterlim:
430 raise RuntimeError("tridiag_eigen: no convergence to an eigenvalue after %d iterations" % iterlim)
432 j += 1
434 # form shift
436 p = d[l]
437 g = (d[l + 1] - p) / (2 * e[l])
438 r = ctx.hypot(g, 1)
440 if g < 0:
441 s = g - r
442 else:
443 s = g + r
445 g = d[m] - p + e[l] / s
447 s, c, p = 1, 1, 0
449 for i in xrange(m - 1, l - 1, -1):
450 f = s * e[i]
451 b = c * e[i]
452 if abs(f) > abs(g): # this here is a slight improvement also used in gaussq.f or acm algorithm 726.
453 c = g / f
454 r = ctx.hypot(c, 1)
455 e[i + 1] = f * r
456 s = 1 / r
457 c = c * s
458 else:
459 s = f / g
460 r = ctx.hypot(s, 1)
461 e[i + 1] = g * r
462 c = 1 / r
463 s = s * c
464 g = d[i + 1] - p
465 r = (d[i] - g) * s + 2 * c * b
466 p = s * r
467 d[i + 1] = g + p
468 g = c * r - b
470 if not isinstance(z, bool):
471 # calculate eigenvectors
472 for w in xrange(z.rows):
473 f = z[w,i+1]
474 z[w,i+1] = s * z[w,i] + c * f
475 z[w,i ] = c * z[w,i] - s * f
477 d[l] = d[l] - p
478 e[l] = g
479 e[m] = 0
481 for ii in xrange(1, n):
482 # sort eigenvalues and eigenvectors (bubble-sort)
483 i = ii - 1
484 k = i
485 p = d[i]
486 for j in xrange(ii, n):
487 if d[j] >= p:
488 continue
489 k = j
490 p = d[k]
491 if k == i:
492 continue
493 d[k] = d[i]
494 d[i] = p
496 if not isinstance(z, bool):
497 for w in xrange(z.rows):
498 p = z[w,i]
499 z[w,i] = z[w,k]
500 z[w,k] = p
502########################################################################################
504@defun
505def eigsy(ctx, A, eigvals_only = False, overwrite_a = False):
506 """
507 This routine solves the (ordinary) eigenvalue problem for a real symmetric
508 square matrix A. Given A, an orthogonal matrix Q is calculated which
509 diagonalizes A:
511 Q' A Q = diag(E) and Q Q' = Q' Q = 1
513 Here diag(E) is a diagonal matrix whose diagonal is E.
514 ' denotes the transpose.
516 The columns of Q are the eigenvectors of A and E contains the eigenvalues:
518 A Q[:,i] = E[i] Q[:,i]
521 input:
523 A: real matrix of format (n,n) which is symmetric
524 (i.e. A=A' or A[i,j]=A[j,i])
526 eigvals_only: if true, calculates only the eigenvalues E.
527 if false, calculates both eigenvectors and eigenvalues.
529 overwrite_a: if true, allows modification of A which may improve
530 performance. if false, A is not modified.
532 output:
534 E: vector of format (n). contains the eigenvalues of A in ascending order.
536 Q: orthogonal matrix of format (n,n). contains the eigenvectors
537 of A as columns.
539 return value:
541 E if eigvals_only is true
542 (E, Q) if eigvals_only is false
544 example:
545 >>> from mpmath import mp
546 >>> A = mp.matrix([[3, 2], [2, 0]])
547 >>> E = mp.eigsy(A, eigvals_only = True)
548 >>> print(E)
549 [-1.0]
550 [ 4.0]
552 >>> A = mp.matrix([[1, 2], [2, 3]])
553 >>> E, Q = mp.eigsy(A)
554 >>> print(mp.chop(A * Q[:,0] - E[0] * Q[:,0]))
555 [0.0]
556 [0.0]
558 see also: eighe, eigh, eig
559 """
561 if not overwrite_a:
562 A = A.copy()
564 d = ctx.zeros(A.rows, 1)
565 e = ctx.zeros(A.rows, 1)
567 if eigvals_only:
568 r_sy_tridiag(ctx, A, d, e, calc_ev = False)
569 tridiag_eigen(ctx, d, e, False)
570 return d
571 else:
572 r_sy_tridiag(ctx, A, d, e, calc_ev = True)
573 tridiag_eigen(ctx, d, e, A)
574 return (d, A)
577@defun
578def eighe(ctx, A, eigvals_only = False, overwrite_a = False):
579 """
580 This routine solves the (ordinary) eigenvalue problem for a complex
581 hermitian square matrix A. Given A, an unitary matrix Q is calculated which
582 diagonalizes A:
584 Q' A Q = diag(E) and Q Q' = Q' Q = 1
586 Here diag(E) a is diagonal matrix whose diagonal is E.
587 ' denotes the hermitian transpose (i.e. ordinary transposition and
588 complex conjugation).
590 The columns of Q are the eigenvectors of A and E contains the eigenvalues:
592 A Q[:,i] = E[i] Q[:,i]
595 input:
597 A: complex matrix of format (n,n) which is hermitian
598 (i.e. A=A' or A[i,j]=conj(A[j,i]))
600 eigvals_only: if true, calculates only the eigenvalues E.
601 if false, calculates both eigenvectors and eigenvalues.
603 overwrite_a: if true, allows modification of A which may improve
604 performance. if false, A is not modified.
606 output:
608 E: vector of format (n). contains the eigenvalues of A in ascending order.
610 Q: unitary matrix of format (n,n). contains the eigenvectors
611 of A as columns.
613 return value:
615 E if eigvals_only is true
616 (E, Q) if eigvals_only is false
618 example:
619 >>> from mpmath import mp
620 >>> A = mp.matrix([[1, -3 - 1j], [-3 + 1j, -2]])
621 >>> E = mp.eighe(A, eigvals_only = True)
622 >>> print(E)
623 [-4.0]
624 [ 3.0]
626 >>> A = mp.matrix([[1, 2 + 5j], [2 - 5j, 3]])
627 >>> E, Q = mp.eighe(A)
628 >>> print(mp.chop(A * Q[:,0] - E[0] * Q[:,0]))
629 [0.0]
630 [0.0]
632 see also: eigsy, eigh, eig
633 """
635 if not overwrite_a:
636 A = A.copy()
638 d = ctx.zeros(A.rows, 1)
639 e = ctx.zeros(A.rows, 1)
640 t = ctx.zeros(A.rows, 1)
642 if eigvals_only:
643 c_he_tridiag_0(ctx, A, d, e, t)
644 tridiag_eigen(ctx, d, e, False)
645 return d
646 else:
647 c_he_tridiag_0(ctx, A, d, e, t)
648 B = ctx.eye(A.rows)
649 tridiag_eigen(ctx, d, e, B)
650 c_he_tridiag_2(ctx, A, t, B)
651 return (d, B)
653@defun
654def eigh(ctx, A, eigvals_only = False, overwrite_a = False):
655 """
656 "eigh" is a unified interface for "eigsy" and "eighe". Depending on
657 whether A is real or complex the appropriate function is called.
659 This routine solves the (ordinary) eigenvalue problem for a real symmetric
660 or complex hermitian square matrix A. Given A, an orthogonal (A real) or
661 unitary (A complex) matrix Q is calculated which diagonalizes A:
663 Q' A Q = diag(E) and Q Q' = Q' Q = 1
665 Here diag(E) a is diagonal matrix whose diagonal is E.
666 ' denotes the hermitian transpose (i.e. ordinary transposition and
667 complex conjugation).
669 The columns of Q are the eigenvectors of A and E contains the eigenvalues:
671 A Q[:,i] = E[i] Q[:,i]
673 input:
675 A: a real or complex square matrix of format (n,n) which is symmetric
676 (i.e. A[i,j]=A[j,i]) or hermitian (i.e. A[i,j]=conj(A[j,i])).
678 eigvals_only: if true, calculates only the eigenvalues E.
679 if false, calculates both eigenvectors and eigenvalues.
681 overwrite_a: if true, allows modification of A which may improve
682 performance. if false, A is not modified.
684 output:
686 E: vector of format (n). contains the eigenvalues of A in ascending order.
688 Q: an orthogonal or unitary matrix of format (n,n). contains the
689 eigenvectors of A as columns.
691 return value:
693 E if eigvals_only is true
694 (E, Q) if eigvals_only is false
696 example:
697 >>> from mpmath import mp
698 >>> A = mp.matrix([[3, 2], [2, 0]])
699 >>> E = mp.eigh(A, eigvals_only = True)
700 >>> print(E)
701 [-1.0]
702 [ 4.0]
704 >>> A = mp.matrix([[1, 2], [2, 3]])
705 >>> E, Q = mp.eigh(A)
706 >>> print(mp.chop(A * Q[:,0] - E[0] * Q[:,0]))
707 [0.0]
708 [0.0]
710 >>> A = mp.matrix([[1, 2 + 5j], [2 - 5j, 3]])
711 >>> E, Q = mp.eigh(A)
712 >>> print(mp.chop(A * Q[:,0] - E[0] * Q[:,0]))
713 [0.0]
714 [0.0]
716 see also: eigsy, eighe, eig
717 """
719 iscomplex = any(type(x) is ctx.mpc for x in A)
721 if iscomplex:
722 return ctx.eighe(A, eigvals_only = eigvals_only, overwrite_a = overwrite_a)
723 else:
724 return ctx.eigsy(A, eigvals_only = eigvals_only, overwrite_a = overwrite_a)
727@defun
728def gauss_quadrature(ctx, n, qtype = "legendre", alpha = 0, beta = 0):
729 """
730 This routine calulates gaussian quadrature rules for different
731 families of orthogonal polynomials. Let (a, b) be an interval,
732 W(x) a positive weight function and n a positive integer.
733 Then the purpose of this routine is to calculate pairs (x_k, w_k)
734 for k=0, 1, 2, ... (n-1) which give
736 int(W(x) * F(x), x = a..b) = sum(w_k * F(x_k),k = 0..(n-1))
738 exact for all polynomials F(x) of degree (strictly) less than 2*n. For all
739 integrable functions F(x) the sum is a (more or less) good approximation to
740 the integral. The x_k are called nodes (which are the zeros of the
741 related orthogonal polynomials) and the w_k are called the weights.
743 parameters
744 n (input) The degree of the quadrature rule, i.e. its number of
745 nodes.
747 qtype (input) The family of orthogonal polynmomials for which to
748 compute the quadrature rule. See the list below.
750 alpha (input) real number, used as parameter for some orthogonal
751 polynomials
753 beta (input) real number, used as parameter for some orthogonal
754 polynomials.
756 return value
758 (X, W) a pair of two real arrays where x_k = X[k] and w_k = W[k].
761 orthogonal polynomials:
763 qtype polynomial
764 ----- ----------
766 "legendre" Legendre polynomials, W(x)=1 on the interval (-1, +1)
767 "legendre01" shifted Legendre polynomials, W(x)=1 on the interval (0, +1)
768 "hermite" Hermite polynomials, W(x)=exp(-x*x) on (-infinity,+infinity)
769 "laguerre" Laguerre polynomials, W(x)=exp(-x) on (0,+infinity)
770 "glaguerre" generalized Laguerre polynomials, W(x)=exp(-x)*x**alpha
771 on (0, +infinity)
772 "chebyshev1" Chebyshev polynomials of the first kind, W(x)=1/sqrt(1-x*x)
773 on (-1, +1)
774 "chebyshev2" Chebyshev polynomials of the second kind, W(x)=sqrt(1-x*x)
775 on (-1, +1)
776 "jacobi" Jacobi polynomials, W(x)=(1-x)**alpha * (1+x)**beta on (-1, +1)
777 with alpha>-1 and beta>-1
779 examples:
780 >>> from mpmath import mp
781 >>> f = lambda x: x**8 + 2 * x**6 - 3 * x**4 + 5 * x**2 - 7
782 >>> X, W = mp.gauss_quadrature(5, "hermite")
783 >>> A = mp.fdot([(f(x), w) for x, w in zip(X, W)])
784 >>> B = mp.sqrt(mp.pi) * 57 / 16
785 >>> C = mp.quad(lambda x: mp.exp(- x * x) * f(x), [-mp.inf, +mp.inf])
786 >>> mp.nprint((mp.chop(A-B, tol = 1e-10), mp.chop(A-C, tol = 1e-10)))
787 (0.0, 0.0)
789 >>> f = lambda x: x**5 - 2 * x**4 + 3 * x**3 - 5 * x**2 + 7 * x - 11
790 >>> X, W = mp.gauss_quadrature(3, "laguerre")
791 >>> A = mp.fdot([(f(x), w) for x, w in zip(X, W)])
792 >>> B = 76
793 >>> C = mp.quad(lambda x: mp.exp(-x) * f(x), [0, +mp.inf])
794 >>> mp.nprint(mp.chop(A-B, tol = 1e-10), mp.chop(A-C, tol = 1e-10))
795 .0
797 # orthogonality of the chebyshev polynomials:
798 >>> f = lambda x: mp.chebyt(3, x) * mp.chebyt(2, x)
799 >>> X, W = mp.gauss_quadrature(3, "chebyshev1")
800 >>> A = mp.fdot([(f(x), w) for x, w in zip(X, W)])
801 >>> print(mp.chop(A, tol = 1e-10))
802 0.0
804 references:
805 - golub and welsch, "calculations of gaussian quadrature rules", mathematics of
806 computation 23, p. 221-230 (1969)
807 - golub, "some modified matrix eigenvalue problems", siam review 15, p. 318-334 (1973)
808 - stroud and secrest, "gaussian quadrature formulas", prentice-hall (1966)
810 See also the routine gaussq.f in netlog.org or ACM Transactions on
811 Mathematical Software algorithm 726.
812 """
814 d = ctx.zeros(n, 1)
815 e = ctx.zeros(n, 1)
816 z = ctx.zeros(1, n)
818 z[0,0] = 1
820 if qtype == "legendre":
821 # legendre on the range -1 +1 , abramowitz, table 25.4, p.916
822 w = 2
823 for i in xrange(n):
824 j = i + 1
825 e[i] = ctx.sqrt(j * j / (4 * j * j - ctx.mpf(1)))
826 elif qtype == "legendre01":
827 # legendre shifted to 0 1 , abramowitz, table 25.8, p.921
828 w = 1
829 for i in xrange(n):
830 d[i] = 1 / ctx.mpf(2)
831 j = i + 1
832 e[i] = ctx.sqrt(j * j / (16 * j * j - ctx.mpf(4)))
833 elif qtype == "hermite":
834 # hermite on the range -inf +inf , abramowitz, table 25.10,p.924
835 w = ctx.sqrt(ctx.pi)
836 for i in xrange(n):
837 j = i + 1
838 e[i] = ctx.sqrt(j / ctx.mpf(2))
839 elif qtype == "laguerre":
840 # laguerre on the range 0 +inf , abramowitz, table 25.9, p. 923
841 w = 1
842 for i in xrange(n):
843 j = i + 1
844 d[i] = 2 * j - 1
845 e[i] = j
846 elif qtype=="chebyshev1":
847 # chebyshev polynimials of the first kind
848 w = ctx.pi
849 for i in xrange(n):
850 e[i] = 1 / ctx.mpf(2)
851 e[0] = ctx.sqrt(1 / ctx.mpf(2))
852 elif qtype == "chebyshev2":
853 # chebyshev polynimials of the second kind
854 w = ctx.pi / 2
855 for i in xrange(n):
856 e[i] = 1 / ctx.mpf(2)
857 elif qtype == "glaguerre":
858 # generalized laguerre on the range 0 +inf
859 w = ctx.gamma(1 + alpha)
860 for i in xrange(n):
861 j = i + 1
862 d[i] = 2 * j - 1 + alpha
863 e[i] = ctx.sqrt(j * (j + alpha))
864 elif qtype == "jacobi":
865 # jacobi polynomials
866 alpha = ctx.mpf(alpha)
867 beta = ctx.mpf(beta)
868 ab = alpha + beta
869 abi = ab + 2
870 w = (2**(ab+1)) * ctx.gamma(alpha + 1) * ctx.gamma(beta + 1) / ctx.gamma(abi)
871 d[0] = (beta - alpha) / abi
872 e[0] = ctx.sqrt(4 * (1 + alpha) * (1 + beta) / ((abi + 1) * (abi * abi)))
873 a2b2 = beta * beta - alpha * alpha
874 for i in xrange(1, n):
875 j = i + 1
876 abi = 2 * j + ab
877 d[i] = a2b2 / ((abi - 2) * abi)
878 e[i] = ctx.sqrt(4 * j * (j + alpha) * (j + beta) * (j + ab) / ((abi * abi - 1) * abi * abi))
879 elif isinstance(qtype, str):
880 raise ValueError("unknown quadrature rule \"%s\"" % qtype)
881 elif not isinstance(qtype, str):
882 w = qtype(d, e)
883 else:
884 assert 0
886 tridiag_eigen(ctx, d, e, z)
888 for i in xrange(len(z)):
889 z[i] *= z[i]
891 z = z.transpose()
892 return (d, w * z)
894##################################################################################################
895##################################################################################################
896##################################################################################################
898def svd_r_raw(ctx, A, V = False, calc_u = False):
899 """
900 This routine computes the singular value decomposition of a matrix A.
901 Given A, two orthogonal matrices U and V are calculated such that
903 A = U S V
905 where S is a suitable shaped matrix whose off-diagonal elements are zero.
906 The diagonal elements of S are the singular values of A, i.e. the
907 squareroots of the eigenvalues of A' A or A A'. Here ' denotes the transpose.
908 Householder bidiagonalization and a variant of the QR algorithm is used.
910 overview of the matrices :
912 A : m*n A gets replaced by U
913 U : m*n U replaces A. If n>m then only the first m*m block of U is
914 non-zero. column-orthogonal: U' U = B
915 here B is a n*n matrix whose first min(m,n) diagonal
916 elements are 1 and all other elements are zero.
917 S : n*n diagonal matrix, only the diagonal elements are stored in
918 the array S. only the first min(m,n) diagonal elements are non-zero.
919 V : n*n orthogonal: V V' = V' V = 1
921 parameters:
922 A (input/output) On input, A contains a real matrix of shape m*n.
923 On output, if calc_u is true A contains the column-orthogonal
924 matrix U; otherwise A is simply used as workspace and thus destroyed.
926 V (input/output) if false, the matrix V is not calculated. otherwise
927 V must be a matrix of shape n*n.
929 calc_u (input) If true, the matrix U is calculated and replaces A.
930 if false, U is not calculated and A is simply destroyed
932 return value:
933 S an array of length n containing the singular values of A sorted by
934 decreasing magnitude. only the first min(m,n) elements are non-zero.
936 This routine is a python translation of the fortran routine svd.f in the
937 software library EISPACK (see netlib.org) which itself is based on the
938 algol procedure svd described in:
939 - num. math. 14, 403-420(1970) by golub and reinsch.
940 - wilkinson/reinsch: handbook for auto. comp., vol ii-linear algebra, 134-151(1971).
942 """
944 m, n = A.rows, A.cols
946 S = ctx.zeros(n, 1)
948 # work is a temporary array of size n
949 work = ctx.zeros(n, 1)
951 g = scale = anorm = 0
952 maxits = 3 * ctx.dps
954 for i in xrange(n): # householder reduction to bidiagonal form
955 work[i] = scale*g
956 g = s = scale = 0
957 if i < m:
958 for k in xrange(i, m):
959 scale += ctx.fabs(A[k,i])
960 if scale != 0:
961 for k in xrange(i, m):
962 A[k,i] /= scale
963 s += A[k,i] * A[k,i]
964 f = A[i,i]
965 g = -ctx.sqrt(s)
966 if f < 0:
967 g = -g
968 h = f * g - s
969 A[i,i] = f - g
970 for j in xrange(i+1, n):
971 s = 0
972 for k in xrange(i, m):
973 s += A[k,i] * A[k,j]
974 f = s / h
975 for k in xrange(i, m):
976 A[k,j] += f * A[k,i]
977 for k in xrange(i,m):
978 A[k,i] *= scale
980 S[i] = scale * g
981 g = s = scale = 0
983 if i < m and i != n - 1:
984 for k in xrange(i+1, n):
985 scale += ctx.fabs(A[i,k])
986 if scale:
987 for k in xrange(i+1, n):
988 A[i,k] /= scale
989 s += A[i,k] * A[i,k]
990 f = A[i,i+1]
991 g = -ctx.sqrt(s)
992 if f < 0:
993 g = -g
994 h = f * g - s
995 A[i,i+1] = f - g
997 for k in xrange(i+1, n):
998 work[k] = A[i,k] / h
1000 for j in xrange(i+1, m):
1001 s = 0
1002 for k in xrange(i+1, n):
1003 s += A[j,k] * A[i,k]
1004 for k in xrange(i+1, n):
1005 A[j,k] += s * work[k]
1007 for k in xrange(i+1, n):
1008 A[i,k] *= scale
1010 anorm = max(anorm, ctx.fabs(S[i]) + ctx.fabs(work[i]))
1012 if not isinstance(V, bool):
1013 for i in xrange(n-2, -1, -1): # accumulation of right hand transformations
1014 V[i+1,i+1] = 1
1016 if work[i+1] != 0:
1017 for j in xrange(i+1, n):
1018 V[i,j] = (A[i,j] / A[i,i+1]) / work[i+1]
1019 for j in xrange(i+1, n):
1020 s = 0
1021 for k in xrange(i+1, n):
1022 s += A[i,k] * V[j,k]
1023 for k in xrange(i+1, n):
1024 V[j,k] += s * V[i,k]
1026 for j in xrange(i+1, n):
1027 V[j,i] = V[i,j] = 0
1029 V[0,0] = 1
1031 if m<n : minnm = m
1032 else : minnm = n
1034 if calc_u:
1035 for i in xrange(minnm-1, -1, -1): # accumulation of left hand transformations
1036 g = S[i]
1037 for j in xrange(i+1, n):
1038 A[i,j] = 0
1039 if g != 0:
1040 g = 1 / g
1041 for j in xrange(i+1, n):
1042 s = 0
1043 for k in xrange(i+1, m):
1044 s += A[k,i] * A[k,j]
1045 f = (s / A[i,i]) * g
1046 for k in xrange(i, m):
1047 A[k,j] += f * A[k,i]
1048 for j in xrange(i, m):
1049 A[j,i] *= g
1050 else:
1051 for j in xrange(i, m):
1052 A[j,i] = 0
1053 A[i,i] += 1
1055 for k in xrange(n - 1, -1, -1):
1056 # diagonalization of the bidiagonal form:
1057 # loop over singular values, and over allowed itations
1059 its = 0
1060 while 1:
1061 its += 1
1062 flag = True
1064 for l in xrange(k, -1, -1):
1065 nm = l-1
1067 if ctx.fabs(work[l]) + anorm == anorm:
1068 flag = False
1069 break
1071 if ctx.fabs(S[nm]) + anorm == anorm:
1072 break
1074 if flag:
1075 c = 0
1076 s = 1
1077 for i in xrange(l, k + 1):
1078 f = s * work[i]
1079 work[i] *= c
1080 if ctx.fabs(f) + anorm == anorm:
1081 break
1082 g = S[i]
1083 h = ctx.hypot(f, g)
1084 S[i] = h
1085 h = 1 / h
1086 c = g * h
1087 s = - f * h
1089 if calc_u:
1090 for j in xrange(m):
1091 y = A[j,nm]
1092 z = A[j,i]
1093 A[j,nm] = y * c + z * s
1094 A[j,i] = z * c - y * s
1096 z = S[k]
1098 if l == k: # convergence
1099 if z < 0: # singular value is made nonnegative
1100 S[k] = -z
1101 if not isinstance(V, bool):
1102 for j in xrange(n):
1103 V[k,j] = -V[k,j]
1104 break
1106 if its >= maxits:
1107 raise RuntimeError("svd: no convergence to an eigenvalue after %d iterations" % its)
1109 x = S[l] # shift from bottom 2 by 2 minor
1110 nm = k-1
1111 y = S[nm]
1112 g = work[nm]
1113 h = work[k]
1114 f = ((y - z) * (y + z) + (g - h) * (g + h))/(2 * h * y)
1115 g = ctx.hypot(f, 1)
1116 if f >= 0: f = ((x - z) * (x + z) + h * ((y / (f + g)) - h)) / x
1117 else: f = ((x - z) * (x + z) + h * ((y / (f - g)) - h)) / x
1119 c = s = 1 # next qt transformation
1121 for j in xrange(l, nm + 1):
1122 g = work[j+1]
1123 y = S[j+1]
1124 h = s * g
1125 g = c * g
1126 z = ctx.hypot(f, h)
1127 work[j] = z
1128 c = f / z
1129 s = h / z
1130 f = x * c + g * s
1131 g = g * c - x * s
1132 h = y * s
1133 y *= c
1134 if not isinstance(V, bool):
1135 for jj in xrange(n):
1136 x = V[j ,jj]
1137 z = V[j+1,jj]
1138 V[j ,jj]= x * c + z * s
1139 V[j+1 ,jj]= z * c - x * s
1140 z = ctx.hypot(f, h)
1141 S[j] = z
1142 if z != 0: # rotation can be arbitray if z=0
1143 z = 1 / z
1144 c = f * z
1145 s = h * z
1146 f = c * g + s * y
1147 x = c * y - s * g
1149 if calc_u:
1150 for jj in xrange(m):
1151 y = A[jj,j ]
1152 z = A[jj,j+1]
1153 A[jj,j ] = y * c + z * s
1154 A[jj,j+1 ] = z * c - y * s
1156 work[l] = 0
1157 work[k] = f
1158 S[k] = x
1160 ##########################
1162 # Sort singular values into decreasing order (bubble-sort)
1164 for i in xrange(n):
1165 imax = i
1166 s = ctx.fabs(S[i]) # s is the current maximal element
1168 for j in xrange(i + 1, n):
1169 c = ctx.fabs(S[j])
1170 if c > s:
1171 s = c
1172 imax = j
1174 if imax != i:
1175 # swap singular values
1177 z = S[i]
1178 S[i] = S[imax]
1179 S[imax] = z
1181 if calc_u:
1182 for j in xrange(m):
1183 z = A[j,i]
1184 A[j,i] = A[j,imax]
1185 A[j,imax] = z
1187 if not isinstance(V, bool):
1188 for j in xrange(n):
1189 z = V[i,j]
1190 V[i,j] = V[imax,j]
1191 V[imax,j] = z
1193 return S
1195#######################
1197def svd_c_raw(ctx, A, V = False, calc_u = False):
1198 """
1199 This routine computes the singular value decomposition of a matrix A.
1200 Given A, two unitary matrices U and V are calculated such that
1202 A = U S V
1204 where S is a suitable shaped matrix whose off-diagonal elements are zero.
1205 The diagonal elements of S are the singular values of A, i.e. the
1206 squareroots of the eigenvalues of A' A or A A'. Here ' denotes the hermitian
1207 transpose (i.e. transposition and conjugation). Householder bidiagonalization
1208 and a variant of the QR algorithm is used.
1210 overview of the matrices :
1212 A : m*n A gets replaced by U
1213 U : m*n U replaces A. If n>m then only the first m*m block of U is
1214 non-zero. column-unitary: U' U = B
1215 here B is a n*n matrix whose first min(m,n) diagonal
1216 elements are 1 and all other elements are zero.
1217 S : n*n diagonal matrix, only the diagonal elements are stored in
1218 the array S. only the first min(m,n) diagonal elements are non-zero.
1219 V : n*n unitary: V V' = V' V = 1
1221 parameters:
1222 A (input/output) On input, A contains a complex matrix of shape m*n.
1223 On output, if calc_u is true A contains the column-unitary
1224 matrix U; otherwise A is simply used as workspace and thus destroyed.
1226 V (input/output) if false, the matrix V is not calculated. otherwise
1227 V must be a matrix of shape n*n.
1229 calc_u (input) If true, the matrix U is calculated and replaces A.
1230 if false, U is not calculated and A is simply destroyed
1232 return value:
1233 S an array of length n containing the singular values of A sorted by
1234 decreasing magnitude. only the first min(m,n) elements are non-zero.
1236 This routine is a python translation of the fortran routine svd.f in the
1237 software library EISPACK (see netlib.org) which itself is based on the
1238 algol procedure svd described in:
1239 - num. math. 14, 403-420(1970) by golub and reinsch.
1240 - wilkinson/reinsch: handbook for auto. comp., vol ii-linear algebra, 134-151(1971).
1242 """
1244 m, n = A.rows, A.cols
1246 S = ctx.zeros(n, 1)
1248 # work is a temporary array of size n
1249 work = ctx.zeros(n, 1)
1250 lbeta = ctx.zeros(n, 1)
1251 rbeta = ctx.zeros(n, 1)
1252 dwork = ctx.zeros(n, 1)
1254 g = scale = anorm = 0
1255 maxits = 3 * ctx.dps
1257 for i in xrange(n): # householder reduction to bidiagonal form
1258 dwork[i] = scale * g # dwork are the side-diagonal elements
1259 g = s = scale = 0
1260 if i < m:
1261 for k in xrange(i, m):
1262 scale += ctx.fabs(ctx.re(A[k,i])) + ctx.fabs(ctx.im(A[k,i]))
1263 if scale != 0:
1264 for k in xrange(i, m):
1265 A[k,i] /= scale
1266 ar = ctx.re(A[k,i])
1267 ai = ctx.im(A[k,i])
1268 s += ar * ar + ai * ai
1269 f = A[i,i]
1270 g = -ctx.sqrt(s)
1271 if ctx.re(f) < 0:
1272 beta = -g - ctx.conj(f)
1273 g = -g
1274 else:
1275 beta = -g + ctx.conj(f)
1276 beta /= ctx.conj(beta)
1277 beta += 1
1278 h = 2 * (ctx.re(f) * g - s)
1279 A[i,i] = f - g
1280 beta /= h
1281 lbeta[i] = (beta / scale) / scale
1282 for j in xrange(i+1, n):
1283 s = 0
1284 for k in xrange(i, m):
1285 s += ctx.conj(A[k,i]) * A[k,j]
1286 f = beta * s
1287 for k in xrange(i, m):
1288 A[k,j] += f * A[k,i]
1289 for k in xrange(i, m):
1290 A[k,i] *= scale
1292 S[i] = scale * g # S are the diagonal elements
1293 g = s = scale = 0
1295 if i < m and i != n - 1:
1296 for k in xrange(i+1, n):
1297 scale += ctx.fabs(ctx.re(A[i,k])) + ctx.fabs(ctx.im(A[i,k]))
1298 if scale:
1299 for k in xrange(i+1, n):
1300 A[i,k] /= scale
1301 ar = ctx.re(A[i,k])
1302 ai = ctx.im(A[i,k])
1303 s += ar * ar + ai * ai
1304 f = A[i,i+1]
1305 g = -ctx.sqrt(s)
1306 if ctx.re(f) < 0:
1307 beta = -g - ctx.conj(f)
1308 g = -g
1309 else:
1310 beta = -g + ctx.conj(f)
1312 beta /= ctx.conj(beta)
1313 beta += 1
1315 h = 2 * (ctx.re(f) * g - s)
1316 A[i,i+1] = f - g
1318 beta /= h
1319 rbeta[i] = (beta / scale) / scale
1321 for k in xrange(i+1, n):
1322 work[k] = A[i, k]
1324 for j in xrange(i+1, m):
1325 s = 0
1326 for k in xrange(i+1, n):
1327 s += ctx.conj(A[i,k]) * A[j,k]
1328 f = s * beta
1329 for k in xrange(i+1,n):
1330 A[j,k] += f * work[k]
1332 for k in xrange(i+1, n):
1333 A[i,k] *= scale
1335 anorm = max(anorm,ctx.fabs(S[i]) + ctx.fabs(dwork[i]))
1337 if not isinstance(V, bool):
1338 for i in xrange(n-2, -1, -1): # accumulation of right hand transformations
1339 V[i+1,i+1] = 1
1341 if dwork[i+1] != 0:
1342 f = ctx.conj(rbeta[i])
1343 for j in xrange(i+1, n):
1344 V[i,j] = A[i,j] * f
1345 for j in xrange(i+1, n):
1346 s = 0
1347 for k in xrange(i+1, n):
1348 s += ctx.conj(A[i,k]) * V[j,k]
1349 for k in xrange(i+1, n):
1350 V[j,k] += s * V[i,k]
1352 for j in xrange(i+1,n):
1353 V[j,i] = V[i,j] = 0
1355 V[0,0] = 1
1357 if m < n : minnm = m
1358 else : minnm = n
1360 if calc_u:
1361 for i in xrange(minnm-1, -1, -1): # accumulation of left hand transformations
1362 g = S[i]
1363 for j in xrange(i+1, n):
1364 A[i,j] = 0
1365 if g != 0:
1366 g = 1 / g
1367 for j in xrange(i+1, n):
1368 s = 0
1369 for k in xrange(i+1, m):
1370 s += ctx.conj(A[k,i]) * A[k,j]
1371 f = s * ctx.conj(lbeta[i])
1372 for k in xrange(i, m):
1373 A[k,j] += f * A[k,i]
1374 for j in xrange(i, m):
1375 A[j,i] *= g
1376 else:
1377 for j in xrange(i, m):
1378 A[j,i] = 0
1379 A[i,i] += 1
1381 for k in xrange(n-1, -1, -1):
1382 # diagonalization of the bidiagonal form:
1383 # loop over singular values, and over allowed itations
1385 its = 0
1386 while 1:
1387 its += 1
1388 flag = True
1390 for l in xrange(k, -1, -1):
1391 nm = l - 1
1393 if ctx.fabs(dwork[l]) + anorm == anorm:
1394 flag = False
1395 break
1397 if ctx.fabs(S[nm]) + anorm == anorm:
1398 break
1400 if flag:
1401 c = 0
1402 s = 1
1403 for i in xrange(l, k+1):
1404 f = s * dwork[i]
1405 dwork[i] *= c
1406 if ctx.fabs(f) + anorm == anorm:
1407 break
1408 g = S[i]
1409 h = ctx.hypot(f, g)
1410 S[i] = h
1411 h = 1 / h
1412 c = g * h
1413 s = -f * h
1415 if calc_u:
1416 for j in xrange(m):
1417 y = A[j,nm]
1418 z = A[j,i]
1419 A[j,nm]= y * c + z * s
1420 A[j,i] = z * c - y * s
1422 z = S[k]
1424 if l == k: # convergence
1425 if z < 0: # singular value is made nonnegative
1426 S[k] = -z
1427 if not isinstance(V, bool):
1428 for j in xrange(n):
1429 V[k,j] = -V[k,j]
1430 break
1432 if its >= maxits:
1433 raise RuntimeError("svd: no convergence to an eigenvalue after %d iterations" % its)
1435 x = S[l] # shift from bottom 2 by 2 minor
1436 nm = k-1
1437 y = S[nm]
1438 g = dwork[nm]
1439 h = dwork[k]
1440 f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2 * h * y)
1441 g = ctx.hypot(f, 1)
1442 if f >=0: f = (( x - z) *( x + z) + h *((y / (f + g)) - h)) / x
1443 else: f = (( x - z) *( x + z) + h *((y / (f - g)) - h)) / x
1445 c = s = 1 # next qt transformation
1447 for j in xrange(l, nm + 1):
1448 g = dwork[j+1]
1449 y = S[j+1]
1450 h = s * g
1451 g = c * g
1452 z = ctx.hypot(f, h)
1453 dwork[j] = z
1454 c = f / z
1455 s = h / z
1456 f = x * c + g * s
1457 g = g * c - x * s
1458 h = y * s
1459 y *= c
1460 if not isinstance(V, bool):
1461 for jj in xrange(n):
1462 x = V[j ,jj]
1463 z = V[j+1,jj]
1464 V[j ,jj]= x * c + z * s
1465 V[j+1,jj ]= z * c - x * s
1466 z = ctx.hypot(f, h)
1467 S[j] = z
1468 if z != 0: # rotation can be arbitray if z=0
1469 z = 1 / z
1470 c = f * z
1471 s = h * z
1472 f = c * g + s * y
1473 x = c * y - s * g
1474 if calc_u:
1475 for jj in xrange(m):
1476 y = A[jj,j ]
1477 z = A[jj,j+1]
1478 A[jj,j ]= y * c + z * s
1479 A[jj,j+1 ]= z * c - y * s
1481 dwork[l] = 0
1482 dwork[k] = f
1483 S[k] = x
1485 ##########################
1487 # Sort singular values into decreasing order (bubble-sort)
1489 for i in xrange(n):
1490 imax = i
1491 s = ctx.fabs(S[i]) # s is the current maximal element
1493 for j in xrange(i + 1, n):
1494 c = ctx.fabs(S[j])
1495 if c > s:
1496 s = c
1497 imax = j
1499 if imax != i:
1500 # swap singular values
1502 z = S[i]
1503 S[i] = S[imax]
1504 S[imax] = z
1506 if calc_u:
1507 for j in xrange(m):
1508 z = A[j,i]
1509 A[j,i] = A[j,imax]
1510 A[j,imax] = z
1512 if not isinstance(V, bool):
1513 for j in xrange(n):
1514 z = V[i,j]
1515 V[i,j] = V[imax,j]
1516 V[imax,j] = z
1518 return S
1520##################################################################################################
1522@defun
1523def svd_r(ctx, A, full_matrices = False, compute_uv = True, overwrite_a = False):
1524 """
1525 This routine computes the singular value decomposition of a matrix A.
1526 Given A, two orthogonal matrices U and V are calculated such that
1528 A = U S V and U' U = 1 and V V' = 1
1530 where S is a suitable shaped matrix whose off-diagonal elements are zero.
1531 Here ' denotes the transpose. The diagonal elements of S are the singular
1532 values of A, i.e. the squareroots of the eigenvalues of A' A or A A'.
1534 input:
1535 A : a real matrix of shape (m, n)
1536 full_matrices : if true, U and V are of shape (m, m) and (n, n).
1537 if false, U and V are of shape (m, min(m, n)) and (min(m, n), n).
1538 compute_uv : if true, U and V are calculated. if false, only S is calculated.
1539 overwrite_a : if true, allows modification of A which may improve
1540 performance. if false, A is not modified.
1542 output:
1543 U : an orthogonal matrix: U' U = 1. if full_matrices is true, U is of
1544 shape (m, m). ortherwise it is of shape (m, min(m, n)).
1546 S : an array of length min(m, n) containing the singular values of A sorted by
1547 decreasing magnitude.
1549 V : an orthogonal matrix: V V' = 1. if full_matrices is true, V is of
1550 shape (n, n). ortherwise it is of shape (min(m, n), n).
1552 return value:
1554 S if compute_uv is false
1555 (U, S, V) if compute_uv is true
1557 overview of the matrices:
1559 full_matrices true:
1560 A : m*n
1561 U : m*m U' U = 1
1562 S as matrix : m*n
1563 V : n*n V V' = 1
1565 full_matrices false:
1566 A : m*n
1567 U : m*min(n,m) U' U = 1
1568 S as matrix : min(m,n)*min(m,n)
1569 V : min(m,n)*n V V' = 1
1571 examples:
1573 >>> from mpmath import mp
1574 >>> A = mp.matrix([[2, -2, -1], [3, 4, -2], [-2, -2, 0]])
1575 >>> S = mp.svd_r(A, compute_uv = False)
1576 >>> print(S)
1577 [6.0]
1578 [3.0]
1579 [1.0]
1581 >>> U, S, V = mp.svd_r(A)
1582 >>> print(mp.chop(A - U * mp.diag(S) * V))
1583 [0.0 0.0 0.0]
1584 [0.0 0.0 0.0]
1585 [0.0 0.0 0.0]
1588 see also: svd, svd_c
1589 """
1591 m, n = A.rows, A.cols
1593 if not compute_uv:
1594 if not overwrite_a:
1595 A = A.copy()
1596 S = svd_r_raw(ctx, A, V = False, calc_u = False)
1597 S = S[:min(m,n)]
1598 return S
1600 if full_matrices and n < m:
1601 V = ctx.zeros(m, m)
1602 A0 = ctx.zeros(m, m)
1603 A0[:,:n] = A
1604 S = svd_r_raw(ctx, A0, V, calc_u = True)
1606 S = S[:n]
1607 V = V[:n,:n]
1609 return (A0, S, V)
1610 else:
1611 if not overwrite_a:
1612 A = A.copy()
1613 V = ctx.zeros(n, n)
1614 S = svd_r_raw(ctx, A, V, calc_u = True)
1616 if n > m:
1617 if full_matrices == False:
1618 V = V[:m,:]
1620 S = S[:m]
1621 A = A[:,:m]
1623 return (A, S, V)
1625##############################
1627@defun
1628def svd_c(ctx, A, full_matrices = False, compute_uv = True, overwrite_a = False):
1629 """
1630 This routine computes the singular value decomposition of a matrix A.
1631 Given A, two unitary matrices U and V are calculated such that
1633 A = U S V and U' U = 1 and V V' = 1
1635 where S is a suitable shaped matrix whose off-diagonal elements are zero.
1636 Here ' denotes the hermitian transpose (i.e. transposition and complex
1637 conjugation). The diagonal elements of S are the singular values of A,
1638 i.e. the squareroots of the eigenvalues of A' A or A A'.
1640 input:
1641 A : a complex matrix of shape (m, n)
1642 full_matrices : if true, U and V are of shape (m, m) and (n, n).
1643 if false, U and V are of shape (m, min(m, n)) and (min(m, n), n).
1644 compute_uv : if true, U and V are calculated. if false, only S is calculated.
1645 overwrite_a : if true, allows modification of A which may improve
1646 performance. if false, A is not modified.
1648 output:
1649 U : an unitary matrix: U' U = 1. if full_matrices is true, U is of
1650 shape (m, m). ortherwise it is of shape (m, min(m, n)).
1652 S : an array of length min(m, n) containing the singular values of A sorted by
1653 decreasing magnitude.
1655 V : an unitary matrix: V V' = 1. if full_matrices is true, V is of
1656 shape (n, n). ortherwise it is of shape (min(m, n), n).
1658 return value:
1660 S if compute_uv is false
1661 (U, S, V) if compute_uv is true
1663 overview of the matrices:
1665 full_matrices true:
1666 A : m*n
1667 U : m*m U' U = 1
1668 S as matrix : m*n
1669 V : n*n V V' = 1
1671 full_matrices false:
1672 A : m*n
1673 U : m*min(n,m) U' U = 1
1674 S as matrix : min(m,n)*min(m,n)
1675 V : min(m,n)*n V V' = 1
1677 example:
1678 >>> from mpmath import mp
1679 >>> A = mp.matrix([[-2j, -1-3j, -2+2j], [2-2j, -1-3j, 1], [-3+1j,-2j,0]])
1680 >>> S = mp.svd_c(A, compute_uv = False)
1681 >>> print(mp.chop(S - mp.matrix([mp.sqrt(34), mp.sqrt(15), mp.sqrt(6)])))
1682 [0.0]
1683 [0.0]
1684 [0.0]
1686 >>> U, S, V = mp.svd_c(A)
1687 >>> print(mp.chop(A - U * mp.diag(S) * V))
1688 [0.0 0.0 0.0]
1689 [0.0 0.0 0.0]
1690 [0.0 0.0 0.0]
1692 see also: svd, svd_r
1693 """
1695 m, n = A.rows, A.cols
1697 if not compute_uv:
1698 if not overwrite_a:
1699 A = A.copy()
1700 S = svd_c_raw(ctx, A, V = False, calc_u = False)
1701 S = S[:min(m,n)]
1702 return S
1704 if full_matrices and n < m:
1705 V = ctx.zeros(m, m)
1706 A0 = ctx.zeros(m, m)
1707 A0[:,:n] = A
1708 S = svd_c_raw(ctx, A0, V, calc_u = True)
1710 S = S[:n]
1711 V = V[:n,:n]
1713 return (A0, S, V)
1714 else:
1715 if not overwrite_a:
1716 A = A.copy()
1717 V = ctx.zeros(n, n)
1718 S = svd_c_raw(ctx, A, V, calc_u = True)
1720 if n > m:
1721 if full_matrices == False:
1722 V = V[:m,:]
1724 S = S[:m]
1725 A = A[:,:m]
1727 return (A, S, V)
1729@defun
1730def svd(ctx, A, full_matrices = False, compute_uv = True, overwrite_a = False):
1731 """
1732 "svd" is a unified interface for "svd_r" and "svd_c". Depending on
1733 whether A is real or complex the appropriate function is called.
1735 This routine computes the singular value decomposition of a matrix A.
1736 Given A, two orthogonal (A real) or unitary (A complex) matrices U and V
1737 are calculated such that
1739 A = U S V and U' U = 1 and V V' = 1
1741 where S is a suitable shaped matrix whose off-diagonal elements are zero.
1742 Here ' denotes the hermitian transpose (i.e. transposition and complex
1743 conjugation). The diagonal elements of S are the singular values of A,
1744 i.e. the squareroots of the eigenvalues of A' A or A A'.
1746 input:
1747 A : a real or complex matrix of shape (m, n)
1748 full_matrices : if true, U and V are of shape (m, m) and (n, n).
1749 if false, U and V are of shape (m, min(m, n)) and (min(m, n), n).
1750 compute_uv : if true, U and V are calculated. if false, only S is calculated.
1751 overwrite_a : if true, allows modification of A which may improve
1752 performance. if false, A is not modified.
1754 output:
1755 U : an orthogonal or unitary matrix: U' U = 1. if full_matrices is true, U is of
1756 shape (m, m). ortherwise it is of shape (m, min(m, n)).
1758 S : an array of length min(m, n) containing the singular values of A sorted by
1759 decreasing magnitude.
1761 V : an orthogonal or unitary matrix: V V' = 1. if full_matrices is true, V is of
1762 shape (n, n). ortherwise it is of shape (min(m, n), n).
1764 return value:
1766 S if compute_uv is false
1767 (U, S, V) if compute_uv is true
1769 overview of the matrices:
1771 full_matrices true:
1772 A : m*n
1773 U : m*m U' U = 1
1774 S as matrix : m*n
1775 V : n*n V V' = 1
1777 full_matrices false:
1778 A : m*n
1779 U : m*min(n,m) U' U = 1
1780 S as matrix : min(m,n)*min(m,n)
1781 V : min(m,n)*n V V' = 1
1783 examples:
1785 >>> from mpmath import mp
1786 >>> A = mp.matrix([[2, -2, -1], [3, 4, -2], [-2, -2, 0]])
1787 >>> S = mp.svd(A, compute_uv = False)
1788 >>> print(S)
1789 [6.0]
1790 [3.0]
1791 [1.0]
1793 >>> U, S, V = mp.svd(A)
1794 >>> print(mp.chop(A - U * mp.diag(S) * V))
1795 [0.0 0.0 0.0]
1796 [0.0 0.0 0.0]
1797 [0.0 0.0 0.0]
1799 see also: svd_r, svd_c
1800 """
1802 iscomplex = any(type(x) is ctx.mpc for x in A)
1804 if iscomplex:
1805 return ctx.svd_c(A, full_matrices = full_matrices, compute_uv = compute_uv, overwrite_a = overwrite_a)
1806 else:
1807 return ctx.svd_r(A, full_matrices = full_matrices, compute_uv = compute_uv, overwrite_a = overwrite_a)