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

1import numpy as np 

2from scipy.linalg import svd 

3 

4 

5__all__ = ['polar'] 

6 

7 

8def polar(a, side="right"): 

9 """ 

10 Compute the polar decomposition. 

11 

12 Returns the factors of the polar decomposition [1]_ `u` and `p` such 

13 that ``a = up`` (if `side` is "right") or ``a = pu`` (if `side` is 

14 "left"), where `p` is positive semidefinite. Depending on the shape 

15 of `a`, either the rows or columns of `u` are orthonormal. When `a` 

16 is a square array, `u` is a square unitary array. When `a` is not 

17 square, the "canonical polar decomposition" [2]_ is computed. 

18 

19 Parameters 

20 ---------- 

21 a : (m, n) array_like 

22 The array to be factored. 

23 side : {'left', 'right'}, optional 

24 Determines whether a right or left polar decomposition is computed. 

25 If `side` is "right", then ``a = up``. If `side` is "left", then 

26 ``a = pu``. The default is "right". 

27 

28 Returns 

29 ------- 

30 u : (m, n) ndarray 

31 If `a` is square, then `u` is unitary. If m > n, then the columns 

32 of `a` are orthonormal, and if m < n, then the rows of `u` are 

33 orthonormal. 

34 p : ndarray 

35 `p` is Hermitian positive semidefinite. If `a` is nonsingular, `p` 

36 is positive definite. The shape of `p` is (n, n) or (m, m), depending 

37 on whether `side` is "right" or "left", respectively. 

38 

39 References 

40 ---------- 

41 .. [1] R. A. Horn and C. R. Johnson, "Matrix Analysis", Cambridge 

42 University Press, 1985. 

43 .. [2] N. J. Higham, "Functions of Matrices: Theory and Computation", 

44 SIAM, 2008. 

45 

46 Examples 

47 -------- 

48 >>> import numpy as np 

49 >>> from scipy.linalg import polar 

50 >>> a = np.array([[1, -1], [2, 4]]) 

51 >>> u, p = polar(a) 

52 >>> u 

53 array([[ 0.85749293, -0.51449576], 

54 [ 0.51449576, 0.85749293]]) 

55 >>> p 

56 array([[ 1.88648444, 1.2004901 ], 

57 [ 1.2004901 , 3.94446746]]) 

58 

59 A non-square example, with m < n: 

60 

61 >>> b = np.array([[0.5, 1, 2], [1.5, 3, 4]]) 

62 >>> u, p = polar(b) 

63 >>> u 

64 array([[-0.21196618, -0.42393237, 0.88054056], 

65 [ 0.39378971, 0.78757942, 0.4739708 ]]) 

66 >>> p 

67 array([[ 0.48470147, 0.96940295, 1.15122648], 

68 [ 0.96940295, 1.9388059 , 2.30245295], 

69 [ 1.15122648, 2.30245295, 3.65696431]]) 

70 >>> u.dot(p) # Verify the decomposition. 

71 array([[ 0.5, 1. , 2. ], 

72 [ 1.5, 3. , 4. ]]) 

73 >>> u.dot(u.T) # The rows of u are orthonormal. 

74 array([[ 1.00000000e+00, -2.07353665e-17], 

75 [ -2.07353665e-17, 1.00000000e+00]]) 

76 

77 Another non-square example, with m > n: 

78 

79 >>> c = b.T 

80 >>> u, p = polar(c) 

81 >>> u 

82 array([[-0.21196618, 0.39378971], 

83 [-0.42393237, 0.78757942], 

84 [ 0.88054056, 0.4739708 ]]) 

85 >>> p 

86 array([[ 1.23116567, 1.93241587], 

87 [ 1.93241587, 4.84930602]]) 

88 >>> u.dot(p) # Verify the decomposition. 

89 array([[ 0.5, 1.5], 

90 [ 1. , 3. ], 

91 [ 2. , 4. ]]) 

92 >>> u.T.dot(u) # The columns of u are orthonormal. 

93 array([[ 1.00000000e+00, -1.26363763e-16], 

94 [ -1.26363763e-16, 1.00000000e+00]]) 

95 

96 """ 

97 if side not in ['right', 'left']: 

98 raise ValueError("`side` must be either 'right' or 'left'") 

99 a = np.asarray(a) 

100 if a.ndim != 2: 

101 raise ValueError("`a` must be a 2-D array.") 

102 

103 w, s, vh = svd(a, full_matrices=False) 

104 u = w.dot(vh) 

105 if side == 'right': 

106 # a = up 

107 p = (vh.T.conj() * s).dot(vh) 

108 else: 

109 # a = pu 

110 p = (w * s).dot(w.T.conj()) 

111 return u, p