Source Code for Module linear_operators.iterative.norms

  1  """ 
  2  Implement norm classes. Available now : 
  3   
  4  - Norm2 
  5  - Huber norm 
  6  - Norm-k 
  7  """ 
  8  import numpy as np 
  9  from ..operators import identity 
 10   
 11  try: 
 12      from scipy.linalg.fblas import dnrm2 
 13  except ImportError: 
 14      pass 
 15   
 16  if 'dnrm2' in locals(): 
 17 -    def norm2(x): 
 18          return dnrm2(x) ** 2 
 19  else: 
 20 -    def norm2(x): 
 21          return np.dot(x.ravel().T, x.ravel()) 
 22   
 23 -def dnorm2(x): 
 24      return 2 * x 
 25   
 26 -def normp(p=2): 
 27      def norm(t): 
 28          return np.sum(np.abs(t) ** p) 
 29      return norm 
 30   
 31 -def dnormp(p=2): 
 32      def norm(t): 
 33          return np.sign(t) * p * (np.abs(t) ** (p - 1)) 
 34      return norm 
 35   
 36 -def hnorm(d=None): 
 37      if d is None: 
 38          return norm2 
 39      else: 
 40          def norm(t): 
 41              return np.sum(huber(t, d)) 
 42          return norm 
 43   
 44 -def dhnorm(d=None): 
 45      if d is None: 
 46          return dnorm2 
 47      else: 
 48          def norm(t): 
 49              return dhuber(t, d) 
 50          return norm 
 51   
 52 -def huber(t, delta=1): 
 53      """Apply the huber function to the vector t, with transition delta""" 
 54      t_out = t.flatten() 
 55      quadratic_index = np.where(np.abs(t_out) < delta) 
 56      linear_index = np.where(np.abs(t_out) >= delta) 
 57      t_out[quadratic_index] = np.abs(t_out[quadratic_index]) ** 2 
 58      t_out[linear_index] = 2 * delta * np.abs(t_out[linear_index]) - delta ** 2 
 59      return np.reshape(t_out, t.shape) 
 60   
 61 -def dhuber(t, delta=1): 
 62      """Apply the derivation of the Huber function to t, transition: delta""" 
 63      t_out = t.flatten() 
 64      quadratic_index = np.where(np.abs(t_out) < delta) 
 65      linear_index_positive = np.where(t_out >= delta) 
 66      linear_index_negative = np.where(t_out <= - delta) 
 67      t_out[quadratic_index] = 2 * t_out[quadratic_index] 
 68      t_out[linear_index_positive] = 2 * delta 
 69      t_out[linear_index_negative] = - 2 * delta 
 70      return np.reshape(t_out, t.shape) 
 71   
 72 -class Norm(object): 
 73      """ 
 74      An abstract class to define norm classes. 
 75      """ 
 76 -    def __call__(self, x): 
 77          return self._call(x) 
 78 -    def diff(self, x): 
 79          return self._diff(x) 
 80   
 81 -class Norm2(Norm): 
 82      """ 
 83      A norm-2 class. Optionally accepts a covariance matrix C. 
 84      If C is given, the norm would be : np.dot(x.T, C * x). 
 85      Otherwise, it would be norm2(x). 
 86   
 87      Parameters 
 88      ---------- 
 89   
 90      C : LinearOperator (None) 
 91          The covariance matrix of the norm. 
 92   
 93      Returns 
 94      ------- 
 95      Returns a Norm2 instance with a __call__ and a diff method. 
 96      """ 
 97 -    def __init__(self, C=None): 
 98          if C is None: 
 99              def call(x): 
100                  return norm2(x) 
101              def diff(x): 
102                  return 2 * x 
103              def hessian(x): 
104                  return 2 * identity(2 *(x.size,)) 
105          else: 
106              def call(x): 
107                  return np.dot(x.T, C * x) 
108              def diff(x): 
109                  return 2 * C * x 
110              def hessian(x): 
111                  return 2 * C 
112          self.C = C 
113          self._call = call 
114          self._diff = diff 
115   
116 -class Huber(Norm): 
117      """ 
118      An Huber norm class. 
119   
120      Parameters 
121      ---------- 
122   
123      delta: float 
124         The Huber parameter of the norm. 
125         if abs(x_i) is below delta, returns x_i ** 2 
126         else returns 2 * delta * x_i - delta ** 2 
127   
128      Returns 
129      ------- 
130      Returns an Huber instance with a __call__ and a diff method. 
131       """ 
132 -    def __init__(self, delta): 
133          self.delta = delta 
134          self._call = hnorm(d=delta) 
135          self._diff = dhnorm(d=delta) 
136   
137 -class Normp(Norm): 
138      """ 
139      An Norm-p class. 
140   
141      Parameters 
142      ---------- 
143   
144      p: float 
145         The power of the norm. 
146         The norm will be np.sum(np.abs(x) ** p) 
147   
148      Returns 
149      ------- 
150      Returns a Normp instance with a __call__ and a diff method. 
151       """ 
152 -    def __init__(self, p): 
153          self.p = p 
154          self._call = normp(p=p) 
155          self._diff = dnormp(p=p) 
156