agd.Metrics.asym_iso
1import numpy as np 2from .base import Base 3from .riemann import Riemann 4from .isotropic import Isotropic 5from .. import AutomaticDifferentiation as ad 6from .. import LinearParallel as lp 7from ..FiniteDifferences import common_field 8 9 10class AsymIso(Base): 11 r""" 12 A special case of AsymQuad metric taking the form 13 $ 14 F(x) = \sqrt{a^2 |x|^2+sign(a) max(0,<w,x>)^2} 15 $ 16 where $a$ is a field of scalars, and w a vector field. 17 Member fields and __init__ arguments : 18 - a : an array of shape (n1,..,nk) 19 - w : an array of shape (vdim,n1,...,nk) 20 """ 21 def __init__(self,a,w): 22 a,w = (ad.asarray(e) for e in (a,w)) 23 self.a,self.w =common_field((a,w),(0,1)) 24 25 def norm(self,v): 26 v,a,w = common_field((ad.asarray(v),self.a,self.w),(1,0,1)) 27 return np.sqrt(a**2*np.sum(v**2,axis=0) + np.sign(a)*np.maximum(lp.dot_VV(w,v),0.)**2) 28 29 def gradient(self,v): 30 v,a,w = common_field((ad.asarray(v),self.a,self.w),(1,0,1)) 31 g = a**2*v + np.sign(a)*np.maximum(0.,lp.dot_VV(w,v))*w 32 return g/np.sqrt(lp.dot_VV(v,g)) 33 34 def dual(self): 35 a2,s,w2 = self.a**2,np.sign(self.a),self._norm2w() 36 r2 = s*(1/a2-1/(a2+s*w2))/w2 37 return AsymIso(-1/self.a,np.sqrt(r2)*self.w) 38 39 @property 40 def vdim(self): return len(self.w) 41 42 @property 43 def shape(self): return self.a.shape 44 45 def _norm2w(self): return np.sum(self.w**2,axis=0) 46 def is_definite(self): 47 return np.where(self.a>0,1.,self.a**2-self._norm2w()) > 0 48 def anisotropy(self): 49 a2,w2 = self.a**2,self._norm2w() 50 return np.sqrt(np.where(self.a>0,1.+w2/a2,1./(1.-w2/a2))) 51 def cost_bound(self): 52 return np.where(self.a>0,np.sqrt(self.a**2+self._norm2w()),-self.a) 53 54 def rotate(self,r): return AsymIso(self.a,lp.dot_AV(r,self.w)) 55 def with_cost(self,cost): return AsymIso(cost*self.a,cost*self.w) 56 def with_costs(self,costs): 57 if len(costs)>1 and not np.allclose(costs[1:],costs[0]): 58 raise ValueError("Costs must be identical along all axes for Metrics.Isotropic. Consider using Metrics.Diagonal class") 59 return self.with_cost(costs[0]) 60 61 def model_HFM(self): 62 return 'AsymIso'+str(self.vdim) 63 64 @classmethod 65 def from_cast(cls,metric): 66 if isinstance(metric,cls): return metric 67 iso = Isotropic.from_cast(metric) 68 w = np.zeros_like(iso.cost,shape=(iso.vdim,)+iso.shape) 69 return cls(iso.cost,w) 70 71 def __iter__(self): 72 yield self.a 73 yield self.w 74 75 def make_proj_dual(self,**kwargs): 76 """kwargs : passed to Riemann.make_proj_dual""" 77 vdim,a = self.vdim,self.a 78 proj_a = Isotropic(np.abs(a)).make_proj_dual() 79 eye = np.eye(vdim,like=a).reshape( (vdim,vdim)+(1,)*a.ndim ) 80 proj_w=Riemann(a**2*eye+np.sign(a)*lp.outer_self(self.w)).make_proj_dual(**kwargs) 81 82 def proj(x): 83 x,w = common_field((x,self.w),depths=(1,1)) 84 x_a = proj_a(x) 85 x_w = proj_w(x) 86 s = lp.dot_VV(w,x_a) > 0 87 return np.where(s[None],x_w,x_a) 88 return proj 89 90 91 92
11class AsymIso(Base): 12 r""" 13 A special case of AsymQuad metric taking the form 14 $ 15 F(x) = \sqrt{a^2 |x|^2+sign(a) max(0,<w,x>)^2} 16 $ 17 where $a$ is a field of scalars, and w a vector field. 18 Member fields and __init__ arguments : 19 - a : an array of shape (n1,..,nk) 20 - w : an array of shape (vdim,n1,...,nk) 21 """ 22 def __init__(self,a,w): 23 a,w = (ad.asarray(e) for e in (a,w)) 24 self.a,self.w =common_field((a,w),(0,1)) 25 26 def norm(self,v): 27 v,a,w = common_field((ad.asarray(v),self.a,self.w),(1,0,1)) 28 return np.sqrt(a**2*np.sum(v**2,axis=0) + np.sign(a)*np.maximum(lp.dot_VV(w,v),0.)**2) 29 30 def gradient(self,v): 31 v,a,w = common_field((ad.asarray(v),self.a,self.w),(1,0,1)) 32 g = a**2*v + np.sign(a)*np.maximum(0.,lp.dot_VV(w,v))*w 33 return g/np.sqrt(lp.dot_VV(v,g)) 34 35 def dual(self): 36 a2,s,w2 = self.a**2,np.sign(self.a),self._norm2w() 37 r2 = s*(1/a2-1/(a2+s*w2))/w2 38 return AsymIso(-1/self.a,np.sqrt(r2)*self.w) 39 40 @property 41 def vdim(self): return len(self.w) 42 43 @property 44 def shape(self): return self.a.shape 45 46 def _norm2w(self): return np.sum(self.w**2,axis=0) 47 def is_definite(self): 48 return np.where(self.a>0,1.,self.a**2-self._norm2w()) > 0 49 def anisotropy(self): 50 a2,w2 = self.a**2,self._norm2w() 51 return np.sqrt(np.where(self.a>0,1.+w2/a2,1./(1.-w2/a2))) 52 def cost_bound(self): 53 return np.where(self.a>0,np.sqrt(self.a**2+self._norm2w()),-self.a) 54 55 def rotate(self,r): return AsymIso(self.a,lp.dot_AV(r,self.w)) 56 def with_cost(self,cost): return AsymIso(cost*self.a,cost*self.w) 57 def with_costs(self,costs): 58 if len(costs)>1 and not np.allclose(costs[1:],costs[0]): 59 raise ValueError("Costs must be identical along all axes for Metrics.Isotropic. Consider using Metrics.Diagonal class") 60 return self.with_cost(costs[0]) 61 62 def model_HFM(self): 63 return 'AsymIso'+str(self.vdim) 64 65 @classmethod 66 def from_cast(cls,metric): 67 if isinstance(metric,cls): return metric 68 iso = Isotropic.from_cast(metric) 69 w = np.zeros_like(iso.cost,shape=(iso.vdim,)+iso.shape) 70 return cls(iso.cost,w) 71 72 def __iter__(self): 73 yield self.a 74 yield self.w 75 76 def make_proj_dual(self,**kwargs): 77 """kwargs : passed to Riemann.make_proj_dual""" 78 vdim,a = self.vdim,self.a 79 proj_a = Isotropic(np.abs(a)).make_proj_dual() 80 eye = np.eye(vdim,like=a).reshape( (vdim,vdim)+(1,)*a.ndim ) 81 proj_w=Riemann(a**2*eye+np.sign(a)*lp.outer_self(self.w)).make_proj_dual(**kwargs) 82 83 def proj(x): 84 x,w = common_field((x,self.w),depths=(1,1)) 85 x_a = proj_a(x) 86 x_w = proj_w(x) 87 s = lp.dot_VV(w,x_a) > 0 88 return np.where(s[None],x_w,x_a) 89 return proj
A special case of AsymQuad metric taking the form
$
F(x) = \sqrt{a^2 |x|^2+sign(a) max(0,
- a : an array of shape (n1,..,nk)
- w : an array of shape (vdim,n1,...,nk)
26 def norm(self,v): 27 v,a,w = common_field((ad.asarray(v),self.a,self.w),(1,0,1)) 28 return np.sqrt(a**2*np.sum(v**2,axis=0) + np.sign(a)*np.maximum(lp.dot_VV(w,v),0.)**2)
Norm or quasi-norm defined by the class, often denoted $N$ in mathematical formulas. Unless incorrect data is provided, this member function obeys, for all vectors $u,v\in R^d$ and all $\alpha \geq 0$
- $N(u+v) \leq N(u)+N(v)$
- $N(\alpha u) = \alpha N(u)$
- $N(u)\geq 0$ with equality iff $u=0$.
Broadcasting will occur depending on the shape of $v$ and of the class data.
30 def gradient(self,v): 31 v,a,w = common_field((ad.asarray(v),self.a,self.w),(1,0,1)) 32 g = a**2*v + np.sign(a)*np.maximum(0.,lp.dot_VV(w,v))*w 33 return g/np.sqrt(lp.dot_VV(v,g))
Gradient of the norm
defined by the class.
35 def dual(self): 36 a2,s,w2 = self.a**2,np.sign(self.a),self._norm2w() 37 r2 = s*(1/a2-1/(a2+s*w2))/w2 38 return AsymIso(-1/self.a,np.sqrt(r2)*self.w)
Dual norm
, mathematically defined by
$N^*(x) = max\{ < x, y> ; N(y)\leq 1 \}$
Dimensions of the underlying domain.
Expected to be the empty tuple, or a tuple of length vdim
.
49 def anisotropy(self): 50 a2,w2 = self.a**2,self._norm2w() 51 return np.sqrt(np.where(self.a>0,1.+w2/a2,1./(1.-w2/a2)))
Anisotropy ratio of the norm
denoted $N$.
Defined as
$$
\max_{|u| = |v| = 1} \frac {N(u)}{N(v)}.
$$
Upper bound on $N(u)$, for any unit vector $u$, where $N$ is the norm
defined by the class.
55 def rotate(self,r): return AsymIso(self.a,lp.dot_AV(r,self.w))
Rotation of the norm, by a given rotation matrix. The new unit ball is the direct image of the previous one.
56 def with_cost(self,cost): return AsymIso(cost*self.a,cost*self.w)
Produces a norm $N'$ obeying $N'(x) = N(cost*x)$.
57 def with_costs(self,costs): 58 if len(costs)>1 and not np.allclose(costs[1:],costs[0]): 59 raise ValueError("Costs must be identical along all axes for Metrics.Isotropic. Consider using Metrics.Diagonal class") 60 return self.with_cost(costs[0])
Produces a norm $N'$ defined by $$ N'(x) = N(costs * x) $$ where the multiplication is elementwise.
65 @classmethod 66 def from_cast(cls,metric): 67 if isinstance(metric,cls): return metric 68 iso = Isotropic.from_cast(metric) 69 w = np.zeros_like(iso.cost,shape=(iso.vdim,)+iso.shape) 70 return cls(iso.cost,w)
Produces a metric by casting another metric of a compatible type.
76 def make_proj_dual(self,**kwargs): 77 """kwargs : passed to Riemann.make_proj_dual""" 78 vdim,a = self.vdim,self.a 79 proj_a = Isotropic(np.abs(a)).make_proj_dual() 80 eye = np.eye(vdim,like=a).reshape( (vdim,vdim)+(1,)*a.ndim ) 81 proj_w=Riemann(a**2*eye+np.sign(a)*lp.outer_self(self.w)).make_proj_dual(**kwargs) 82 83 def proj(x): 84 x,w = common_field((x,self.w),depths=(1,1)) 85 x_a = proj_a(x) 86 x_w = proj_w(x) 87 s = lp.dot_VV(w,x_a) > 0 88 return np.where(s[None],x_w,x_a) 89 return proj
kwargs : passed to Riemann.make_proj_dual