Coverage for suppy\feasibility\_hyperplanes\_kaczmarz_extrapolations.py: 86%

1from typing import List 

2import numpy as np 

3import numpy.typing as npt 

4from scipy import sparse 

5 

6try: 

7 import cupy as cp 

8 

9 NO_GPU = False 

10 

11except ImportError: 

12 NO_GPU = True 

13 cp = np 

14 

15 

16from suppy.feasibility._hyperplanes._kaczmarz_algorithms import SimultaneousKaczmarzMethod 

17 

18 

19class ExtrapolatedLandweberHyperplane(SimultaneousKaczmarzMethod): 

20 """ 

21 Extrapolated Landweber algorithm for solving linear equalities of the 

22 form Ax = b. 

23 

24 Parameters 

25 ---------- 

26 A : npt.NDArray or sparse.sparray 

27 Matrix for linear systems 

28 b : npt.NDArray 

29 Bound for linear systems 

30 relaxation : float, optional 

31 The outer relaxation parameter, by default 1.0. 

32 weights : None or List[float], optional 

33 The weights for the constraints, by default None. 

34 proximity_flag : bool, optional 

35 Flag to indicate if proximity calculations should be performed, by default True. 

36 

37 References 

38 ---------- 

39 - [1] https://doi.org/10.1007/s11075-025-02163-0 

40 - [2] https://doi.org/10.1007/978-3-642-30901-4 

41 """ 

42 

43 def __init__( 

44 self, 

45 A: npt.NDArray | sparse.sparray, 

46 b: npt.NDArray, 

47 relaxation: float = 1.0, 

48 weights: None | List[float] | npt.NDArray = None, 

49 proximity_flag: bool = True, 

50 ): 

51 super().__init__(A, b, 1.0, relaxation, weights, proximity_flag) 

52 self.a_i = self.A.row_norm(2, 2) 

53 self.weight_norm = self.weights / self.a_i 

54 self.sigmas = [] 

55 

56 def _project(self, x: npt.NDArray) -> np.ndarray: 

57 xp = cp if self._use_gpu else np 

58 p = self.map(x) 

59 res = self.b - p 

60 res_idx = res != 0 

61 if not xp.any(res_idx): 

62 self.sigmas.append(0) 

63 return x 

64 t = self.weight_norm * res 

65 t_2 = t @ self.A 

66 sig = (res @ t) / (t_2 @ t_2) 

67 self.sigmas.append(sig) 

68 x += sig * t_2 

69 

70 return x 

71 

72 

73class AdaptiveStepLandweberHyperplane(SimultaneousKaczmarzMethod): 

74 """ 

75 Landweber algorithm with adaptive step size for solving linear 

76 equalities 

77 of the form Ax = b. 

78 

79 Parameters 

80 ---------- 

81 A : npt.NDArray or sparse.sparray 

82 Matrix for linear systems 

83 b : npt.NDArray 

84 Bound for linear systems 

85 relaxation : float, optional 

86 The outer relaxation parameter, by default 1.0. 

87 weights : None or List[float], optional 

88 The weights for the constraints, by default None. 

89 proximity_flag : bool, optional 

90 Flag to indicate if proximity calculations should be performed, by default True. 

91 

92 References 

93 ---------- 

94 - [1] https://doi.org/10.1515/jiip-2015-0082 

95 """ 

96 

97 def __init__( 

98 self, 

99 A: npt.NDArray | sparse.sparray, 

100 b: npt.NDArray, 

101 relaxation: float = 1.0, 

102 weights: None | List[float] | npt.NDArray = None, 

103 proximity_flag: bool = True, 

104 ): 

105 super().__init__(A, b, 1.0, relaxation, weights, proximity_flag) 

106 self.sigmas = [] 

107 

108 def _project(self, x: npt.NDArray) -> np.ndarray: 

109 xp = cp if self._use_gpu else np 

110 p = self.map(x) 

111 res = self.b - p 

112 res_idx = res != 0 

113 if not xp.any(res_idx): 

114 self.sigmas.append(0) 

115 return x 

116 u_t = (self.weights * res) @ self.A 

117 Au_t = self.A @ u_t 

118 sig = (u_t @ u_t) / (Au_t @ (self.weights * Au_t)) 

119 self.sigmas.append(sig) 

120 x += sig * u_t 

121 

122 return x