agd.Metrics.asym_quad

View Source

  1# Copyright 2020 Jean-Marie Mirebeau, University Paris-Sud, CNRS, University Paris-Saclay
  2# Distributed WITHOUT ANY WARRANTY. Licensed under the Apache License, Version 2.0, see http://www.apache.org/licenses/LICENSE-2.0
  3
  4import numpy as np
  5from .base import Base
  6from .riemann import Riemann
  7from .rander import Rander
  8from . import misc
  9from .. import LinearParallel as lp
 10from .. import AutomaticDifferentiation as ad
 11from ..FiniteDifferences import common_field
 12
 13class AsymQuad(Base):
 14	r"""
 15	An asymmetric quadratic norm takes the form 
 16	$
 17	F(x) = \sqrt{< x, m x > + max(0,< w, x >)^2},
 18	$
 19	where $m$ is a given symmetric positive definite tensor, 
 20	and w is a given vector.
 21
 22	Member fields and __init__ arguments : 
 23	- m : an array of shape (vdim,vdim,n1,..,nk) where vdim is the ambient space dimension.
 24	The array must be symmetric, a.k.a m[i,j] = m[j,i] for all 0 <= i < j < vdim.
 25	- w : an array of shape (vdim,n1,...,nk)
 26	"""
 27
 28	def __init__(self,m,w):
 29		m,w = (ad.asarray(e) for e in (m,w))
 30		self.m,self.w =common_field((m,w),(2,1))
 31
 32	def norm(self,v):
 33		v,m,w = common_field((ad.asarray(v),self.m,self.w),(1,2,1))
 34		return np.sqrt(lp.dot_VAV(v,m,v) + np.maximum(lp.dot_VV(w,v),0.)**2)
 35
 36	def gradient(self,v):
 37		v,m,w = common_field((ad.asarray(v),self.m,self.w),(1,2,1))
 38		a = lp.dot_AV(m,v) + w*np.maximum(0.,lp.dot_VV(w,v))
 39		return a/np.sqrt(lp.dot_VV(v,a))
 40
 41
 42	def dual(self):
 43		M = lp.inverse(self.m+lp.outer_self(self.w))
 44		wInv = lp.solve_AV(self.m,self.w)
 45		W = -wInv/np.sqrt(1.+lp.dot_VV(self.w,wInv))
 46		return AsymQuad(M,W)
 47
 48	@property
 49	def vdim(self): return len(self.m)
 50
 51	@property
 52	def shape(self): return self.m.shape[2:]	
 53
 54	def is_definite(self):
 55		return Riemann(self.m).is_definite()
 56
 57	def anisotropy(self):
 58		eMax = Riemann(self.m+lp.outer_self(self.w)).eigvals().max(axis=0)
 59		eMin = Riemann(self.m).eigvals().min(axis=0)
 60		return np.sqrt(eMax/eMin)
 61
 62	def cost_bound(self):
 63		return Riemann(self.m + lp.outer_self(self.w)).cost_bound()
 64
 65	def inv_transform(self,a):
 66		rander = Rander(self.m,self.w).inv_transform(a)
 67		return AsymQuad(rander.m,rander.w)
 68	def with_costs(self,costs):
 69		rander = Rander(self.m,self.w).with_costs(costs)
 70		return AsymQuad(rander.m,rander.w)
 71
 72	def flatten(self,solve_w=False):
 73		m,w = self.m,self.w
 74		if solve_w: w = lp.solve_AV(m,w)
 75		return Rander(m,w).flatten()
 76
 77	@classmethod
 78	def expand(cls,arr):
 79		rd = Rander.expand(arr)
 80		return cls(rd.m,rd.w)
 81
 82	def model_HFM(self):
 83		return "AsymmetricQuadratic"+str(self.vdim)
 84
 85
 86	@classmethod
 87	def needle(cls,u,cost_forward,cost_orthogonal,cost_reverse=None):
 88		"""
 89		Defines a needle like metric
 90		- u : reference direction. Denote U = u/|u|, and V the orthogonal unit vector.
 91		- cost_forward =  norm(U)
 92		- cost_orthogonal = norm(V) = norm(-V)
 93		- cost_reverse = norm(-U). (Defaults to cost_orthogonal)
 94		"""
 95		if cost_reverse is None: cost_reverse = cost_orthogonal
 96		cost_parallel = np.minimum(cost_forward,cost_reverse)
 97		riem,_u = Riemann.needle(u,cost_parallel,cost_orthogonal,ret_u=True)
 98		cost_diff = cost_forward**2-cost_reverse**2
 99		return cls(riem.m,np.sign(cost_diff)*np.sqrt(np.abs(cost_diff))*_u)
100
101	@classmethod
102	def from_cast(cls,metric): 
103		if isinstance(metric,cls):	return metric
104		if metric.model_HFM().startswith('AsymIso'):
105			vdim,a,s,w = metric.vdim,metric.a,np.sign(metric.a),metric.w
106			eye = np.eye(vdim,like=a).reshape( (vdim,vdim)+(1,)*a.ndim )
107			return cls(a**2*eye-(s<0)*lp.outer_self(w),s*w)
108		riemann = Riemann.from_cast(metric)
109		return cls(riemann.m,(0.,)*riemann.vdim)
110
111	def __iter__(self):
112		yield self.m
113		yield self.w
114
115	def make_proj_dual(self,**kwargs):
116		"""kwargs : passed to Riemann.make_proj_dual"""
117		proj_m = Riemann(self.m).make_proj_dual(**kwargs)
118		proj_w = Riemann(self.m + lp.outer_self(self.w)).make_proj_dual(**kwargs)
119		def proj(x):
120			x,w = common_field((x,self.w),depths=(1,1))
121			x_m = proj_m(x)
122			x_w = proj_w(x)
123			s = lp.dot_VV(w,x_m) > 0
124			return np.where(s[None],x_w,x_m)
125		return proj

class AsymQuad(agd.Metrics.base.Base): View Source

 14class AsymQuad(Base):
 15	r"""
 16	An asymmetric quadratic norm takes the form 
 17	$
 18	F(x) = \sqrt{< x, m x > + max(0,< w, x >)^2},
 19	$
 20	where $m$ is a given symmetric positive definite tensor, 
 21	and w is a given vector.
 22
 23	Member fields and __init__ arguments : 
 24	- m : an array of shape (vdim,vdim,n1,..,nk) where vdim is the ambient space dimension.
 25	The array must be symmetric, a.k.a m[i,j] = m[j,i] for all 0 <= i < j < vdim.
 26	- w : an array of shape (vdim,n1,...,nk)
 27	"""
 28
 29	def __init__(self,m,w):
 30		m,w = (ad.asarray(e) for e in (m,w))
 31		self.m,self.w =common_field((m,w),(2,1))
 32
 33	def norm(self,v):
 34		v,m,w = common_field((ad.asarray(v),self.m,self.w),(1,2,1))
 35		return np.sqrt(lp.dot_VAV(v,m,v) + np.maximum(lp.dot_VV(w,v),0.)**2)
 36
 37	def gradient(self,v):
 38		v,m,w = common_field((ad.asarray(v),self.m,self.w),(1,2,1))
 39		a = lp.dot_AV(m,v) + w*np.maximum(0.,lp.dot_VV(w,v))
 40		return a/np.sqrt(lp.dot_VV(v,a))
 41
 42
 43	def dual(self):
 44		M = lp.inverse(self.m+lp.outer_self(self.w))
 45		wInv = lp.solve_AV(self.m,self.w)
 46		W = -wInv/np.sqrt(1.+lp.dot_VV(self.w,wInv))
 47		return AsymQuad(M,W)
 48
 49	@property
 50	def vdim(self): return len(self.m)
 51
 52	@property
 53	def shape(self): return self.m.shape[2:]	
 54
 55	def is_definite(self):
 56		return Riemann(self.m).is_definite()
 57
 58	def anisotropy(self):
 59		eMax = Riemann(self.m+lp.outer_self(self.w)).eigvals().max(axis=0)
 60		eMin = Riemann(self.m).eigvals().min(axis=0)
 61		return np.sqrt(eMax/eMin)
 62
 63	def cost_bound(self):
 64		return Riemann(self.m + lp.outer_self(self.w)).cost_bound()
 65
 66	def inv_transform(self,a):
 67		rander = Rander(self.m,self.w).inv_transform(a)
 68		return AsymQuad(rander.m,rander.w)
 69	def with_costs(self,costs):
 70		rander = Rander(self.m,self.w).with_costs(costs)
 71		return AsymQuad(rander.m,rander.w)
 72
 73	def flatten(self,solve_w=False):
 74		m,w = self.m,self.w
 75		if solve_w: w = lp.solve_AV(m,w)
 76		return Rander(m,w).flatten()
 77
 78	@classmethod
 79	def expand(cls,arr):
 80		rd = Rander.expand(arr)
 81		return cls(rd.m,rd.w)
 82
 83	def model_HFM(self):
 84		return "AsymmetricQuadratic"+str(self.vdim)
 85
 86
 87	@classmethod
 88	def needle(cls,u,cost_forward,cost_orthogonal,cost_reverse=None):
 89		"""
 90		Defines a needle like metric
 91		- u : reference direction. Denote U = u/|u|, and V the orthogonal unit vector.
 92		- cost_forward =  norm(U)
 93		- cost_orthogonal = norm(V) = norm(-V)
 94		- cost_reverse = norm(-U). (Defaults to cost_orthogonal)
 95		"""
 96		if cost_reverse is None: cost_reverse = cost_orthogonal
 97		cost_parallel = np.minimum(cost_forward,cost_reverse)
 98		riem,_u = Riemann.needle(u,cost_parallel,cost_orthogonal,ret_u=True)
 99		cost_diff = cost_forward**2-cost_reverse**2
100		return cls(riem.m,np.sign(cost_diff)*np.sqrt(np.abs(cost_diff))*_u)
101
102	@classmethod
103	def from_cast(cls,metric): 
104		if isinstance(metric,cls):	return metric
105		if metric.model_HFM().startswith('AsymIso'):
106			vdim,a,s,w = metric.vdim,metric.a,np.sign(metric.a),metric.w
107			eye = np.eye(vdim,like=a).reshape( (vdim,vdim)+(1,)*a.ndim )
108			return cls(a**2*eye-(s<0)*lp.outer_self(w),s*w)
109		riemann = Riemann.from_cast(metric)
110		return cls(riemann.m,(0.,)*riemann.vdim)
111
112	def __iter__(self):
113		yield self.m
114		yield self.w
115
116	def make_proj_dual(self,**kwargs):
117		"""kwargs : passed to Riemann.make_proj_dual"""
118		proj_m = Riemann(self.m).make_proj_dual(**kwargs)
119		proj_w = Riemann(self.m + lp.outer_self(self.w)).make_proj_dual(**kwargs)
120		def proj(x):
121			x,w = common_field((x,self.w),depths=(1,1))
122			x_m = proj_m(x)
123			x_w = proj_w(x)
124			s = lp.dot_VV(w,x_m) > 0
125			return np.where(s[None],x_w,x_m)
126		return proj

An asymmetric quadratic norm takes the form $ F(x) = \sqrt{< x, m x > + max(0,< w, x >)^2}, $ where $m$ is a given symmetric positive definite tensor, and w is a given vector.

Member fields and __init__ arguments :

m : an array of shape (vdim,vdim,n1,..,nk) where vdim is the ambient space dimension. The array must be symmetric, a.k.a m[i,j] = m[j,i] for all 0 <= i < j < vdim.
w : an array of shape (vdim,n1,...,nk)

AsymQuad(m, w) View Source

29	def __init__(self,m,w):
30		m,w = (ad.asarray(e) for e in (m,w))
31		self.m,self.w =common_field((m,w),(2,1))

def norm(self, v): View Source

33	def norm(self,v):
34		v,m,w = common_field((ad.asarray(v),self.m,self.w),(1,2,1))
35		return np.sqrt(lp.dot_VAV(v,m,v) + 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.

def gradient(self, v): View Source

37	def gradient(self,v):
38		v,m,w = common_field((ad.asarray(v),self.m,self.w),(1,2,1))
39		a = lp.dot_AV(m,v) + w*np.maximum(0.,lp.dot_VV(w,v))
40		return a/np.sqrt(lp.dot_VV(v,a))

Gradient of the norm defined by the class.

def dual(self): View Source

43	def dual(self):
44		M = lp.inverse(self.m+lp.outer_self(self.w))
45		wInv = lp.solve_AV(self.m,self.w)
46		W = -wInv/np.sqrt(1.+lp.dot_VV(self.w,wInv))
47		return AsymQuad(M,W)

Dual norm, mathematically defined by $N^*(x) = max\{ < x, y> ; N(y)\leq 1 \}$

vdim View Source

49	@property
50	def vdim(self): return len(self.m)

The ambient vector space dimension, often denoted $d$ in mathematical formulas.

shape View Source

52	@property
53	def shape(self): return self.m.shape[2:]

Dimensions of the underlying domain. Expected to be the empty tuple, or a tuple of length vdim.

def is_definite(self): View Source

55	def is_definite(self):
56		return Riemann(self.m).is_definite()

Attempts to check wether the data defines a mathematically valid norm.

def anisotropy(self): View Source

58	def anisotropy(self):
59		eMax = Riemann(self.m+lp.outer_self(self.w)).eigvals().max(axis=0)
60		eMin = Riemann(self.m).eigvals().min(axis=0)
61		return np.sqrt(eMax/eMin)

Anisotropy ratio of the norm denoted $N$. Defined as $$ \max_{|u| = |v| = 1} \frac {N(u)}{N(v)}. $$

def cost_bound(self): View Source

63	def cost_bound(self):
64		return Riemann(self.m + lp.outer_self(self.w)).cost_bound()

Upper bound on $N(u)$, for any unit vector $u$, where $N$ is the norm defined by the class.

def inv_transform(self, a): View Source

66	def inv_transform(self,a):
67		rander = Rander(self.m,self.w).inv_transform(a)
68		return AsymQuad(rander.m,rander.w)

Affine transformation of the norm. The new unit ball is the inverse image of the previous one.

def with_costs(self, costs): View Source

69	def with_costs(self,costs):
70		rander = Rander(self.m,self.w).with_costs(costs)
71		return AsymQuad(rander.m,rander.w)

Produces a norm $N'$ defined by $$ N'(x) = N(costs * x) $$ where the multiplication is elementwise.

def flatten(self, solve_w=False): View Source

73	def flatten(self,solve_w=False):
74		m,w = self.m,self.w
75		if solve_w: w = lp.solve_AV(m,w)
76		return Rander(m,w).flatten()

Flattens and concatenate the member fields into a single array.

@classmethod

def expand(cls, arr): View Source

78	@classmethod
79	def expand(cls,arr):
80		rd = Rander.expand(arr)
81		return cls(rd.m,rd.w)

Inverse of the flatten member function. Turns a suitable array into a metric.

def model_HFM(self): View Source

83	def model_HFM(self):
84		return "AsymmetricQuadratic"+str(self.vdim)

The name of the 'model' for parameter, as input to the HFM library.

@classmethod

def needle(cls, u, cost_forward, cost_orthogonal, cost_reverse=None): View Source

 87	@classmethod
 88	def needle(cls,u,cost_forward,cost_orthogonal,cost_reverse=None):
 89		"""
 90		Defines a needle like metric
 91		- u : reference direction. Denote U = u/|u|, and V the orthogonal unit vector.
 92		- cost_forward =  norm(U)
 93		- cost_orthogonal = norm(V) = norm(-V)
 94		- cost_reverse = norm(-U). (Defaults to cost_orthogonal)
 95		"""
 96		if cost_reverse is None: cost_reverse = cost_orthogonal
 97		cost_parallel = np.minimum(cost_forward,cost_reverse)
 98		riem,_u = Riemann.needle(u,cost_parallel,cost_orthogonal,ret_u=True)
 99		cost_diff = cost_forward**2-cost_reverse**2
100		return cls(riem.m,np.sign(cost_diff)*np.sqrt(np.abs(cost_diff))*_u)

Defines a needle like metric

u : reference direction. Denote U = u/|u|, and V the orthogonal unit vector.
cost_forward = norm(U)
cost_orthogonal = norm(V) = norm(-V)
cost_reverse = norm(-U). (Defaults to cost_orthogonal)

@classmethod

def from_cast(cls, metric): View Source

102	@classmethod
103	def from_cast(cls,metric): 
104		if isinstance(metric,cls):	return metric
105		if metric.model_HFM().startswith('AsymIso'):
106			vdim,a,s,w = metric.vdim,metric.a,np.sign(metric.a),metric.w
107			eye = np.eye(vdim,like=a).reshape( (vdim,vdim)+(1,)*a.ndim )
108			return cls(a**2*eye-(s<0)*lp.outer_self(w),s*w)
109		riemann = Riemann.from_cast(metric)
110		return cls(riemann.m,(0.,)*riemann.vdim)

Produces a metric by casting another metric of a compatible type.

def make_proj_dual(self, **kwargs): View Source

116	def make_proj_dual(self,**kwargs):
117		"""kwargs : passed to Riemann.make_proj_dual"""
118		proj_m = Riemann(self.m).make_proj_dual(**kwargs)
119		proj_w = Riemann(self.m + lp.outer_self(self.w)).make_proj_dual(**kwargs)
120		def proj(x):
121			x,w = common_field((x,self.w),depths=(1,1))
122			x_m = proj_m(x)
123			x_w = proj_w(x)
124			s = lp.dot_VV(w,x_m) > 0
125			return np.where(s[None],x_w,x_m)
126		return proj

kwargs : passed to Riemann.make_proj_dual

Inherited Members

agd.Metrics.base.Base: disassociate; anisotropy_bound; cos_asym; cos; angle; transform; rotate; rotate_by; with_speeds; with_cost; with_speed; to_HFM; from_HFM; from_generator; array_float_caster; norm2; gradient2; set_interpolation; at; make_proj; make_prox