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

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"""
2Linear Algebra solvers and other helpers
3"""
4import numpy as np
6__all__ = ["logdet_symm", "stationary_solve", "transf_constraints",
7 "matrix_sqrt"]
10def logdet_symm(m, check_symm=False):
11 """
12 Return log(det(m)) asserting positive definiteness of m.
14 Parameters
15 ----------
16 m : array_like
17 2d array that is positive-definite (and symmetric)
19 Returns
20 -------
21 logdet : float
22 The log-determinant of m.
23 """
24 from scipy import linalg
25 if check_symm:
26 if not np.all(m == m.T): # would be nice to short-circuit check
27 raise ValueError("m is not symmetric.")
28 c, _ = linalg.cho_factor(m, lower=True)
29 return 2*np.sum(np.log(c.diagonal()))
32def stationary_solve(r, b):
33 """
34 Solve a linear system for a Toeplitz correlation matrix.
36 A Toeplitz correlation matrix represents the covariance of a
37 stationary series with unit variance.
39 Parameters
40 ----------
41 r : array_like
42 A vector describing the coefficient matrix. r[0] is the first
43 band next to the diagonal, r[1] is the second band, etc.
44 b : array_like
45 The right-hand side for which we are solving, i.e. we solve
46 Tx = b and return b, where T is the Toeplitz coefficient matrix.
48 Returns
49 -------
50 The solution to the linear system.
51 """
53 db = r[0:1]
55 dim = b.ndim
56 if b.ndim == 1:
57 b = b[:, None]
58 x = b[0:1, :]
60 for j in range(1, len(b)):
61 rf = r[0:j][::-1]
62 a = (b[j, :] - np.dot(rf, x)) / (1 - np.dot(rf, db[::-1]))
63 z = x - np.outer(db[::-1], a)
64 x = np.concatenate((z, a[None, :]), axis=0)
66 if j == len(b) - 1:
67 break
69 rn = r[j]
70 a = (rn - np.dot(rf, db)) / (1 - np.dot(rf, db[::-1]))
71 z = db - a*db[::-1]
72 db = np.concatenate((z, np.r_[a]))
74 if dim == 1:
75 x = x[:, 0]
77 return x
80def transf_constraints(constraints):
81 """use QR to get transformation matrix to impose constraint
83 Parameters
84 ----------
85 constraints : ndarray, 2-D
86 restriction matrix with one constraints in rows
88 Returns
89 -------
90 transf : ndarray
91 transformation matrix to reparameterize so that constraint is
92 imposed
94 Notes
95 -----
96 This is currently and internal helper function for GAM.
97 API not stable and will most likely change.
99 The code for this function was taken from patsy spline handling, and
100 corresponds to the reparameterization used by Wood in R's mgcv package.
102 See Also
103 --------
104 statsmodels.base._constraints.TransformRestriction : class to impose
105 constraints by reparameterization used by `_fit_constrained`.
106 """
108 from scipy import linalg
110 m = constraints.shape[0]
111 q, _ = linalg.qr(np.transpose(constraints))
112 transf = q[:, m:]
113 return transf
116def matrix_sqrt(mat, inverse=False, full=False, nullspace=False,
117 threshold=1e-15):
118 """matrix square root for symmetric matrices
120 Usage is for decomposing a covariance function S into a square root R
121 such that
123 R' R = S if inverse is False, or
124 R' R = pinv(S) if inverse is True
126 Parameters
127 ----------
128 mat : array_like, 2-d square
129 symmetric square matrix for which square root or inverse square
130 root is computed.
131 There is no checking for whether the matrix is symmetric.
132 A warning is issued if some singular values are negative, i.e.
133 below the negative of the threshold.
134 inverse : bool
135 If False (default), then the matrix square root is returned.
136 If inverse is True, then the matrix square root of the inverse
137 matrix is returned.
138 full : bool
139 If full is False (default, then the square root has reduce number
140 of rows if the matrix is singular, i.e. has singular values below
141 the threshold.
142 nullspace: bool
143 If nullspace is true, then the matrix square root of the null space
144 of the matrix is returned.
145 threshold : float
146 Singular values below the threshold are dropped.
148 Returns
149 -------
150 msqrt : ndarray
151 matrix square root or square root of inverse matrix.
152 """
153 # see also scipy.linalg null_space
154 u, s, v = np.linalg.svd(mat)
155 if np.any(s < -threshold):
156 import warnings
157 warnings.warn('some singular values are negative')
159 if not nullspace:
160 mask = s > threshold
161 s[s < threshold] = 0
162 else:
163 mask = s < threshold
164 s[s > threshold] = 0
166 sqrt_s = np.sqrt(s[mask])
167 if inverse:
168 sqrt_s = 1 / np.sqrt(s[mask])
170 if full:
171 b = np.dot(u[:, mask], np.dot(np.diag(sqrt_s), v[mask]))
172 else:
173 b = np.dot(np.diag(sqrt_s), v[mask])
174 return b