Coverage for src/driada/utils/matrix.py: 16.67%

1from numpy import linalg as la 

2import numpy as np 

3 

4 

5def nearestPD(A): 

6 """Find the nearest positive-definite matrix to input 

7 

8 A Python/Numpy port of John D'Errico's `nearestSPD` MATLAB code [1], which 

9 credits [2]. 

10 

11 [1] https://www.mathworks.com/matlabcentral/fileexchange/42885-nearestspd 

12 

13 [2] N.J. Higham, "Computing a nearest symmetric positive semidefinite 

14 matrix" (1988): https://doi.org/10.1016/0024-3795(88)90223-6 

15 """ 

16 

17 B = (A + A.T) / 2 

18 _, s, V = la.svd(B) 

19 

20 H = np.dot(V.T, np.dot(np.diag(s), V)) 

21 

22 A2 = (B + H) / 2 

23 

24 A3 = (A2 + A2.T) / 2 

25 

26 if isPD(A3): 

27 return A3 

28 

29 spacing = np.spacing(la.norm(A)) 

30 # The above is different from [1]. It appears that MATLAB's `chol` Cholesky 

31 # decomposition will accept matrixes with exactly 0-eigenvalue, whereas 

32 # Numpy's will not. So where [1] uses `eps(mineig)` (where `eps` is Matlab 

33 # for `np.spacing`), we use the above definition. CAVEAT: our `spacing` 

34 # will be much larger than [1]'s `eps(mineig)`, since `mineig` is usually on 

35 # the order of 1e-16, and `eps(1e-16)` is on the order of 1e-34, whereas 

36 # `spacing` will, for Gaussian random matrixes of small dimension, be on 

37 # the order of 1e-16. 

38 I = np.eye(A.shape[0]) 

39 k = 1 

40 while not isPD(A3): 

41 mineig = np.min(np.real(la.eigvals(A3))) 

42 A3 += I * (-mineig * k**2 + spacing) 

43 k += 1 

44 

45 return A3 

46 

47 

48def isPD(B): 

49 """Returns true when input is positive-definite, via Cholesky""" 

50 try: 

51 _ = la.cholesky(B) 

52 return True 

53 except la.LinAlgError: 

54 return False