Hide keyboard shortcuts

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# -*- coding: utf-8 -*- 

2""" 

3This file contains analytic implementations of rotation methods. 

4""" 

5 

6import numpy as np 

7import scipy as sp 

8 

9 

10def target_rotation(A, H, full_rank=False): 

11 r""" 

12 Analytically performs orthogonal rotations towards a target matrix, 

13 i.e., we minimize: 

14 

15 .. math:: 

16 \phi(L) =\frac{1}{2}\|AT-H\|^2. 

17 

18 where :math:`T` is an orthogonal matrix. This problem is also known as 

19 an orthogonal Procrustes problem. 

20 

21 Under the assumption that :math:`A^*H` has full rank, the analytical 

22 solution :math:`T` is given by: 

23 

24 .. math:: 

25 T = (A^*HH^*A)^{-\frac{1}{2}}A^*H, 

26 

27 see Green (1952). In other cases the solution is given by :math:`T = UV`, 

28 where :math:`U` and :math:`V` result from the singular value decomposition 

29 of :math:`A^*H`: 

30 

31 .. math:: 

32 A^*H = U\Sigma V, 

33 

34 see Schonemann (1966). 

35 

36 Parameters 

37 ---------- 

38 A : numpy matrix (default None) 

39 non rotated factors 

40 H : numpy matrix 

41 target matrix 

42 full_rank : bool (default FAlse) 

43 if set to true full rank is assumed 

44 

45 Returns 

46 ------- 

47 The matrix :math:`T`. 

48 

49 References 

50 ---------- 

51 [1] Green (1952, Psychometrika) - The orthogonal approximation of an 

52 oblique structure in factor analysis 

53 

54 [2] Schonemann (1966) - A generalized solution of the orthogonal 

55 procrustes problem 

56 

57 [3] Gower, Dijksterhuis (2004) - Procrustes problems 

58 """ 

59 ATH = A.T.dot(H) 

60 if full_rank or np.linalg.matrix_rank(ATH) == A.shape[1]: 

61 T = sp.linalg.fractional_matrix_power(ATH.dot(ATH.T), -1/2).dot(ATH) 

62 else: 

63 U, D, V = np.linalg.svd(ATH, full_matrices=False) 

64 T = U.dot(V) 

65 return T 

66 

67 

68def procrustes(A, H): 

69 r""" 

70 Analytically solves the following Procrustes problem: 

71 

72 .. math:: 

73 \phi(L) =\frac{1}{2}\|AT-H\|^2. 

74 

75 (With no further conditions on :math:`H`) 

76 

77 Under the assumption that :math:`A^*H` has full rank, the analytical 

78 solution :math:`T` is given by: 

79 

80 .. math:: 

81 T = (A^*HH^*A)^{-\frac{1}{2}}A^*H, 

82 

83 see Navarra, Simoncini (2010). 

84 

85 Parameters 

86 ---------- 

87 A : numpy matrix 

88 non rotated factors 

89 H : numpy matrix 

90 target matrix 

91 full_rank : bool (default False) 

92 if set to true full rank is assumed 

93 

94 Returns 

95 ------- 

96 The matrix :math:`T`. 

97 

98 References 

99 ---------- 

100 [1] Navarra, Simoncini (2010) - A guide to empirical orthogonal functions 

101 for climate data analysis 

102 """ 

103 return np.linalg.inv(A.T.dot(A)).dot(A.T).dot(H) 

104 

105 

106def promax(A, k=2): 

107 r""" 

108 Performs promax rotation of the matrix :math:`A`. 

109 

110 This method was not very clear to me from the literature, this 

111 implementation is as I understand it should work. 

112 

113 Promax rotation is performed in the following steps: 

114 

115 * Determine varimax rotated patterns :math:`V`. 

116 

117 * Construct a rotation target matrix :math:`|V_{ij}|^k/V_{ij}` 

118 

119 * Perform procrustes rotation towards the target to obtain T 

120 

121 * Determine the patterns 

122 

123 First, varimax rotation a target matrix :math:`H` is determined with 

124 orthogonal varimax rotation. 

125 Then, oblique target rotation is performed towards the target. 

126 

127 Parameters 

128 ---------- 

129 A : numpy matrix 

130 non rotated factors 

131 k : float 

132 parameter, should be positive 

133 

134 References 

135 ---------- 

136 [1] Browne (2001) - An overview of analytic rotation in exploratory 

137 factor analysis 

138 

139 [2] Navarra, Simoncini (2010) - A guide to empirical orthogonal functions 

140 for climate data analysis 

141 """ 

142 assert k > 0 

143 # define rotation target using varimax rotation: 

144 from ._wrappers import rotate_factors 

145 V, T = rotate_factors(A, 'varimax') 

146 H = np.abs(V)**k/V 

147 # solve procrustes problem 

148 S = procrustes(A, H) # np.linalg.inv(A.T.dot(A)).dot(A.T).dot(H); 

149 # normalize 

150 d = np.sqrt(np.diag(np.linalg.inv(S.T.dot(S)))) 

151 D = np.diag(d) 

152 T = np.linalg.inv(S.dot(D)).T 

153 return A.dot(T), T