# Coverage for src/mlstatpy/ml/matrices.py: 100%

## 124 statements

, created at 2022-06-13 20:42 +0200

1# -*- coding: utf-8 -*-

2"""

3@file

5"""

6import warnings

7import numpy

8import numpy.linalg

9from scipy.linalg.lapack import dtrtri # pylint: disable=E0611

12def gram_schmidt(mat, change=False):

13 """

14 Applies the Gram–Schmidt process

15 <https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process>_.

16 Due to performance, every row is considered as a vector.

18 @param mat matrix

19 @param change returns the matrix to change the basis

20 @return new matrix or (new matrix, change matrix)

22 The function assumes the matrix *mat* is

23 horizontal: it has more columns than rows.

25 .. note::

26 The implementation could be improved

27 by directly using :epkg:BLAS function.

29 .. runpython::

30 :showcode:

32 import numpy

33 from mlstatpy.ml.matrices import gram_schmidt

35 X = numpy.array([[1., 2., 3., 4.],

36 [5., 6., 6., 6.],

37 [5., 6., 7., 8.]])

38 T, P = gram_schmidt(X, change=True)

39 print(T)

40 print(P)

41 """

42 if len(mat.shape) != 2:

43 raise ValueError("mat must be a matrix.") # pragma: no cover

44 if mat.shape[1] < mat.shape[0]:

45 raise RuntimeError( # pragma: no cover

46 "The function only works if the number of rows is less "

47 "than the number of columns.")

48 if change:

49 base = numpy.identity(mat.shape[0])

50 # The following code is equivalent to:

51 # res = numpy.empty(mat.shape)

52 # for i in range(0, mat.shape[0]):

53 # res[i, :] = mat[i, :]

54 # for j in range(0, i):

55 # d = numpy.dot(res[j, :], mat[i, :])

56 # res[i, :] -= res[j, :] * d

57 # if change:

58 # base[i, :] -= base[j, :] * d

59 # d = numpy.dot(res[i, :], res[i, :])

60 # if d > 0:

61 # d **= 0.5

62 # res[i, :] /= d

63 # if change:

64 # base[i, :] /= d

65 # But it is faster to write it this way:

66 res = numpy.empty(mat.shape)

67 for i in range(0, mat.shape[0]):

68 res[i, :] = mat[i, :]

69 if i > 0:

70 d = numpy.dot(res[:i, :], mat[i, :])

71 m = numpy.multiply(res[:i, :], d.reshape((len(d), 1)))

72 m = numpy.sum(m, axis=0)

73 res[i, :] -= m

74 if change:

75 m = numpy.multiply(base[:i, :], d.reshape((len(d), 1)))

76 m = numpy.sum(m, axis=0)

77 base[i, :] -= m

78 d = numpy.dot(res[i, :], res[i, :])

79 if d > 0:

80 d **= 0.5

81 res[i, :] /= d

82 if change:

83 base[i, :] /= d

84 return (res, base) if change else res

87def linear_regression(X, y, algo=None):

88 """

89 Solves the linear regression problem,

90 find :math:\\beta which minimizes

91 :math:\\norme{y - X\\beta}, based on the algorithm

92 :ref:Arbre de décision optimisé pour les régressions linéaires

93 <algo_decision_tree_mselin>.

95 @param X features

96 @param y targets

97 @param algo None to use the standard algorithm

98 :math:\\beta = (X'X)^{-1} X'y,

99 'gram', 'qr'

100 @return beta

102 .. runpython::

103 :showcode:

105 import numpy

106 from mlstatpy.ml.matrices import linear_regression

108 X = numpy.array([[1., 2., 3., 4.],

109 [5., 6., 6., 6.],

110 [5., 6., 7., 8.]]).T

111 y = numpy.array([0.1, 0.2, 0.19, 0.29])

112 beta = linear_regression(X, y, algo="gram")

113 print(beta)

115 algo=None computes :math:\\beta = (X'X)^{-1} X'y.

116 algo='qr' uses a QR <https://docs.scipy.org/doc/numpy/reference

117 /generated/numpy.linalg.qr.html>_ decomposition and calls function

118 dtrtri <https://docs.scipy.org/doc/scipy/reference/generated/scipy.

119 linalg.lapack.dtrtri.html>_ to invert an upper triangular matrix.

120 algo='gram' uses :func:gram_schmidt

121 <mlstatpy.ml.matrices.gram_schmidt> and then computes

122 the solution of the linear regression (see above for a link

123 to the algorithm).

124 """

125 if len(y.shape) != 1:

126 warnings.warn( # pragma: no cover

127 "This function is not tested for a multidimensional linear regression.")

128 if algo is None:

129 inv = numpy.linalg.inv(X.T @ X)

130 return inv @ (X.T @ y)

131 if algo == "gram":

132 T, P = gram_schmidt(X.T, change=True)

133 # T = P X

134 return (y.T @ T.T @ P).ravel()

135 if algo == "qr":

136 Q, R = numpy.linalg.qr(X, "full")

137 Ri = dtrtri(R)[0]

138 gamma = (y.T @ Q).ravel()

139 return (gamma @ Ri.T).ravel()

140 raise ValueError( # pragma: no cover

141 "Unknwown algo='{}'.".format(algo))

144def norm2(X):

145 """

146 Computes the square norm for all rows of a

147 matrix.

148 """

149 res = numpy.empty(X.shape[1])

150 for i in range(X.shape[1]):

151 res[i] = numpy.dot(X[i], X[i])

152 return res

155def streaming_gram_schmidt_update(Xk, Pk):

156 """

157 Updates matrix :math:P_k to produce :math:P_{k+1}

158 which is the matrix *P* in algorithm

159 :ref:Streaming Linear Regression

160 <algo_reg_lin_gram_schmidt_streaming>.

161 The function modifies the matrix *Pk*

162 given as an input.

164 @param Xk kth row

165 @param Pk matrix *P* at iteration *k-1*

166 """

167 tki = Pk @ Xk

168 idi = numpy.identity(Pk.shape[0])

170 for i in range(0, Pk.shape[0]):

171 val = tki[i]

173 if i > 0:

174 # for j in range(0, i):

175 # d = tki[j] * val

176 # tki[i] -= tki[j] * d

177 # Pk[i, :] -= Pk[j, :] * d

178 # idi[i, :] -= idi[j, :] * d

180 dv = tki[:i] * val

181 tki[i] -= numpy.dot(dv, tki[:i])

182 dv = dv.reshape((i, 1))

183 Pk[i, :] -= numpy.multiply(Pk[:i, :], dv).sum(axis=0)

184 idi[i, :] -= numpy.multiply(idi[:i, :], dv).sum(axis=0)

186 d = numpy.square(idi[i, :]).sum() # pylint: disable=E1101

187 d = tki[i] ** 2 + d

188 if d > 0:

189 d **= 0.5

190 d = 1. / d

191 tki[i] *= d

192 Pk[i, :] *= d

193 idi[i, :] *= d

196def streaming_gram_schmidt(mat, start=None):

197 """

198 Solves the linear regression problem,

199 find :math:\\beta which minimizes

200 :math:\\norme{y - X\\beta}, based on

201 algorithm :ref:Streaming Gram-Schmidt

202 <algo_reg_lin_gram_schmidt_streaming>.

204 @param mat matrix

205 @param start first row to start iteration, X.shape[1] by default

206 @return iterator on

208 The function assumes the matrix *mat* is

209 horizontal: it has more columns than rows.

211 .. runpython::

212 :showcode:

214 import numpy

215 from mlstatpy.ml.matrices import streaming_gram_schmidt

217 X = numpy.array([[1, 0.5, 10., 5., -2.],

218 [0, 0.4, 20, 4., 2.],

219 [0, 0.7, 20, 4., 2.]], dtype=float).T

221 for i, p in enumerate(streaming_gram_schmidt(X.T)):

222 print("iteration", i, "\\n", p)

223 t = X[:i+3] @ p.T

224 print(t.T @ t)

225 """

226 if len(mat.shape) != 2:

227 raise ValueError("mat must be a matrix.") # pragma: no cover

228 if mat.shape[1] < mat.shape[0]:

229 raise RuntimeError("The function only works if the number of rows is less "

230 "than the number of columns.")

231 if start is None:

232 start = mat.shape[0]

233 mats = mat[:, :start]

234 _, Pk = gram_schmidt(mats, change=True)

235 yield Pk

237 k = start

238 while k < mat.shape[1]:

239 streaming_gram_schmidt_update(mat[:, k], Pk)

240 yield Pk

241 k += 1

244def streaming_linear_regression_update(Xk, yk, XkXk, bk):

245 """

246 Updates coefficients :math:\\beta_k to produce :math:\\beta_{k+1}

247 in :ref:l-piecewise-linear-regression.

248 The function modifies the matrix *Pk*

249 given as an input.

251 @param Xk kth row

252 @param yk kth target

253 @param XkXk matrix :math:X_{1..k}'X_{1..k}', updated by the function

254 @param bk current coefficient (updated by the function)

255 """

256 Xk = Xk.reshape((1, XkXk.shape[0]))

257 xxk = Xk.T @ Xk

258 XkXk += xxk

259 err = Xk.T * (yk - Xk @ bk)

260 bk[:] += (numpy.linalg.inv(XkXk) @ err).flatten()

263def streaming_linear_regression(mat, y, start=None):

264 """

265 Streaming algorithm to solve a linear regression.

266 See :ref:l-piecewise-linear-regression.

268 @param mat features

269 @param y expected target

270 @return iterator on coefficients

272 .. runpython::

273 :showcode:

275 import numpy

276 from mlstatpy.ml.matrices import streaming_linear_regression, linear_regression

278 X = numpy.array([[1, 0.5, 10., 5., -2.],

279 [0, 0.4, 20, 4., 2.],

280 [0, 0.7, 20, 4., 3.]], dtype=float).T

281 y = numpy.array([1., 0.3, 10, 5.1, -3.])

283 for i, bk in enumerate(streaming_linear_regression(X, y)):

284 bk0 = linear_regression(X[:i+3], y[:i+3])

285 print("iteration", i, bk, bk0)

286 """

287 if len(mat.shape) != 2:

288 raise ValueError("mat must be a matrix.") # pragma: no cover

289 if mat.shape[0] < mat.shape[1]:

290 raise RuntimeError("The function only works if the number of rows is more "

291 "than the number of columns.")

292 if len(y.shape) != 1:

293 warnings.warn( # pragma: no cover

294 "This function is not tested for a multidimensional linear regression.")

295 if start is None:

296 start = mat.shape[1]

298 Xk = mat[:start]

299 XkXk = Xk.T @ Xk

300 bk = numpy.linalg.inv(XkXk) @ (Xk.T @ y[:start])

301 yield bk

303 k = start

304 while k < mat.shape[0]:

305 streaming_linear_regression_update(mat[k], y[k], XkXk, bk)

306 yield bk

307 k += 1

310def streaming_linear_regression_gram_schmidt_update(Xk, yk, Xkyk, Pk, bk):

311 """

312 Updates coefficients :math:\\beta_k to produce :math:\\beta_{k+1}

313 in :ref:Streaming Linear Regression

314 <l-piecewise-linear-regression-gram_schmidt>.

315 The function modifies the matrix *Pk*

316 given as an input.

318 @param Xk kth row

319 @param yk kth target

320 @param Xkyk matrix :math:X_{1..k}' y_{1..k}' (updated by the function)

321 @param Pk Gram-Schmidt matrix produced by the streaming algorithm

322 (updated by the function)

323 @param bk current coefficient (updated by the function)

324 """

325 Xk = Xk.T

326 streaming_gram_schmidt_update(Xk, Pk)

327 Xkyk += (Xk * yk).reshape(Xkyk.shape)

328 bk[:] = Pk @ Xkyk @ Pk

331def streaming_linear_regression_gram_schmidt(mat, y, start=None):

332 """

333 Streaming algorithm to solve a linear regression with

334 Gram-Schmidt algorithm.

335 See :ref:l-piecewise-linear-regression-gram_schmidt.

337 @param mat features

338 @param y expected target

339 @return iterator on coefficients

341 .. runpython::

342 :showcode:

344 import numpy

345 from mlstatpy.ml.matrices import streaming_linear_regression, linear_regression

347 X = numpy.array([[1, 0.5, 10., 5., -2.],

348 [0, 0.4, 20, 4., 2.],

349 [0, 0.7, 20, 4., 3.]], dtype=float).T

350 y = numpy.array([1., 0.3, 10, 5.1, -3.])

352 for i, bk in enumerate(streaming_linear_regression(X, y)):

353 bk0 = linear_regression(X[:i+3], y[:i+3])

354 print("iteration", i, bk, bk0)

355 """

356 if len(mat.shape) != 2:

357 raise ValueError("mat must be a matrix.") # pragma: no cover

358 if mat.shape[0] < mat.shape[1]:

359 raise RuntimeError("The function only works if the number of rows is more "

360 "than the number of columns.")

361 if len(y.shape) != 1:

362 warnings.warn( # pragma: no cover

363 "This function is not tested for a multidimensional linear regression.")

364 if start is None:

365 start = mat.shape[1]

367 Xk = mat[:start]

368 xyk = Xk.T @ y[:start]

369 _, Pk = gram_schmidt(Xk.T, change=True)

370 bk = Pk @ xyk @ Pk

371 yield bk

373 k = start

374 while k < mat.shape[0]:

375 streaming_linear_regression_gram_schmidt_update(

376 mat[k], y[k], xyk, Pk, bk)

377 yield bk

378 k += 1