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

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

""" 

This module implements the split-Bregman algorithm described in 

T. Goldstein and S. Osher, SIAM J. Imaging Sci. 2, 323 (2009); 

doi:10.1137/080725891 

""" 

 

import numpy as np 

from scipy.optimize import minimize 

from typing import Any, Dict 

 

 

def fit_split_bregman(A: np.ndarray, 

y: np.ndarray, 

mu: float = 1e-3, 

lmbda: float = 100, 

n_iters: int = 1000, 

tol: float = 1e-6) -> Dict[str, Any]: 

""" 

Determines the solution :math:`\\boldsymbol{x}` to the linear 

problem :math:`\\boldsymbol{A}\\boldsymbol{x}=\\boldsymbol{y}` using 

the split-Bregman algorithm described in T. Goldstein and S. Osher, 

SIAM J. Imaging Sci. 2, 323 (2009); doi:10.1137/080725891. 

The thus obtained parameters are returned in the form of a 

dictionary with a key named `parameters` 

 

Parameters 

---------- 

A 

fit matrix 

y 

target array 

mu 

sparseness parameter 

lmbda 

weight of additional L2-norm in split-Bregman 

n_iters 

maximal number of split-Bregman iterations 

tol 

convergence criterion iterative minimization 

""" 

 

def _shrink(y: np.ndarray, alpha: float) -> np.ndarray: 

""" 

Shrinkage operator as defined in Eq. (11) of the paper by Nelson 

et al., Phys. Rev. B 87, 035125 (2013); doi:10.1103/PhysRevB.87.035125. 

""" 

return np.sign(y) * np.maximum(np.abs(y) - alpha, 0.0) 

 

n_cols = A.shape[1] 

d = np.zeros(n_cols) 

b = np.zeros(n_cols) 

x = np.zeros(n_cols) 

 

old_norm = 0.0 

 

# Precompute for speed. 

AtA = np.dot(A.conj().transpose(), A) 

ftA = np.dot(y.conj().transpose(), A) 

ii = 0 

60 ↛ 78line 60 didn't jump to line 78, because the loop on line 60 didn't complete for i in range(n_iters): 

args = (A, y, mu, lmbda, d, b, AtA, ftA) 

res = minimize(_objective_function, x, args, method='BFGS', 

options={'disp': False}, 

jac=_objective_function_derivative) 

x = res.x 

 

d = _shrink(mu*x + b, 1.0/lmbda) 

b = b + mu*x - d 

 

new_norm = np.linalg.norm(x) 

ii = ii + 1 

 

if abs(new_norm-old_norm) < tol: 

break 

 

old_norm = new_norm 

 

fit_results = {'parameters': x} 

return fit_results 

 

 

def _objective_function(x: np.ndarray, A: np.ndarray, y: np.ndarray, 

mu: float, lmbda: float, d: np.ndarray, b: np.ndarray, 

AtA: np.ndarray, ftA: np.ndarray) -> np.ndarray: 

""" 

Returns the objective function to be minimized. 

 

Parameters 

----------- 

X 

fit matrix 

y 

target array 

mu 

the parameter that adjusts sparseness. 

lmbda 

Split Bregman parameter 

d 

same notation as Nelson et al. paper 

b 

same notation as Nelson et al. paper 

AtA 

sensing matrix transpose times sensing matrix. 

ftA 

np.dot(y.conj().transpose(), A) 

""" 

 

error_vector = np.dot(A, x) - y 

 

obj_function = 0.5*np.vdot(error_vector, error_vector) 

 

112 ↛ 113line 112 didn't jump to line 113, because the condition on line 112 was never true if obj_function.imag > 0.0: 

raise RuntimeError( 

'Objective function contains non-zero imaginary part.)') 

 

sparseness_correction = d - b - mu*x 

obj_function += 0.5*lmbda * \ 

np.vdot(sparseness_correction, sparseness_correction) 

 

120 ↛ 121line 120 didn't jump to line 121, because the condition on line 120 was never true if obj_function.imag > 0.0: 

raise RuntimeError( 

'Objective function contains non-zero imaginary part.)') 

 

return obj_function 

 

 

def _objective_function_derivative(x: np.ndarray, 

A: np.ndarray, 

y: np.ndarray, 

mu: float, 

lmbda: float, 

d: np.ndarray, 

b: np.ndarray, 

AtA: np.ndarray, 

ftA: np.ndarray) -> np.ndarray: 

""" 

Returns the derivative of the objective function. 

 

Parameters 

----------- 

X 

fit matrix 

y 

target array 

mu 

the parameter that adjusts sparseness. 

lmbda 

Split Bregman parameter 

d 

same notation as Nelson, Hart paper 

b 

same notation as Nelson, Hart paper 

AtA 

sensing matrix transpose times sensing matrix. 

ftA 

np.dot(y.conj().transpose(), A) 

""" 

ret = np.squeeze(np.dot(x[np.newaxis, :], AtA) - 

ftA - lmbda*mu*(d - mu * x - b)) 

return ret