agd.Metrics.diagonal

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 .isotropic import Isotropic
 7from .. import AutomaticDifferentiation as ad
 8from .. import FiniteDifferences as fd
 9
10class Diagonal(Base):
11	r"""
12	A Diagonal norm takes the form 
13	$$
14	F(x) = \sqrt{\sum_{0\leq i < d} c_i x_i^2 },
15	$$
16	where $(c_i)_{0\leq i < d}$, are given positive scalars
17
18	Member fields and __init__ arguments : 
19	- costs : the array of $(c_i)_{0\leq i < d}$ positive numbers. Required shape : $(d,n_1,..,n_k)$ where $d=$`vdim` is the ambient space dimension.
20	"""
21
22	def __init__(self,costs):
23		self.costs = ad.asarray(costs)
24
25	@classmethod
26	def from_speed(cls,speeds): 
27		"""
28		Produces a metric whose costs equal 1/speeds
29		"""
30		return cls(1./speeds)
31
32	def dual(self): return self.from_speed(self.costs)
33	def with_costs(self,costs): 
34		self_costs,costs = fd.common_field((self.costs,costs),(1,1))
35		return Diagonal(self_costs*costs)
36
37	def norm(self,v):
38		costs,v = fd.common_field((self.costs,v),depths=(1,1))
39		return ad.Optimization.norm(costs*v,ord=2,axis=0)
40
41	def is_definite(self): return np.all(self.costs>0.,axis=0)
42	def anisotropy(self): return np.max(self.costs,axis=0)/np.min(self.costs,axis=0)
43	def cost_bound(self): return np.max(self.costs,axis=0)
44
45	@property
46	def vdim(self): return len(self.costs)
47	@property
48	def shape(self): return self.costs.shape[1:]
49
50	def flatten(self):      return self.costs
51	@classmethod
52	def expand(cls,arr):    return cls(arr)
53
54	def model_HFM(self):
55		return "Diagonal"+str(self.vdim)
56
57	@classmethod
58	def from_cast(cls,metric):
59		if isinstance(metric,cls): return metric
60		iso = Isotropic.from_cast(metric)
61		shape = (iso.vdim,) + iso.shape
62		return cls(np.broadcast_to(iso.cost,shape))
63
64	def __iter__(self):
65		yield self.costs
66
67	def make_proj_dual(self,newton_maxiter=7):
68		"""
69		Returns the orthogonal projection onto the unit ball.
70		- newton_maxiter : number of iterations for the inner subproblem
71		"""
72		# Implementation based on Benamou, J.-D., Carlier, G. & Hatchi, R. A numerical solution to 
73		# Monge’s problem with a Finsler distance as cost. ESAIM: M2AN 52, 2133–2148 (2018).
74		d_ = self.costs
75		def proj(x):
76			d,x = fd.common_field((d_,x),depths=(1,1))
77			x2=x**2
78			norm2 = np.sum(x2/d,axis=0)
79			# Solve, via a Newton method, the equation f(β)=0 where
80			# f(β) = -1 + sum_i xi^2/(1+di*β)^2 .
81			# By convexity of f, no damping or fancy stopping criterion is needed
82			β = np.min(d,axis=0)*(np.sqrt(norm2)-1) # Initial guess, exact in isotropic case.
83			orig = β<=0.; β[orig]=0. # Those are projected onto the origin
84			lc2 = np.where(orig[None],1.,d*x2)
85			for i in range(newton_maxiter):
86				ilb = 1./(d+β)
87				a = lc2*ilb**2
88				val = np.sum(a,axis=0)-1.
89				val[orig]=0
90				der = -2*np.sum(a*ilb,axis=0)
91				β -= val/der
92			return x*d/(d+β)
93		return proj

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

11class Diagonal(Base):
12	r"""
13	A Diagonal norm takes the form 
14	$$
15	F(x) = \sqrt{\sum_{0\leq i < d} c_i x_i^2 },
16	$$
17	where $(c_i)_{0\leq i < d}$, are given positive scalars
18
19	Member fields and __init__ arguments : 
20	- costs : the array of $(c_i)_{0\leq i < d}$ positive numbers. Required shape : $(d,n_1,..,n_k)$ where $d=$`vdim` is the ambient space dimension.
21	"""
22
23	def __init__(self,costs):
24		self.costs = ad.asarray(costs)
25
26	@classmethod
27	def from_speed(cls,speeds): 
28		"""
29		Produces a metric whose costs equal 1/speeds
30		"""
31		return cls(1./speeds)
32
33	def dual(self): return self.from_speed(self.costs)
34	def with_costs(self,costs): 
35		self_costs,costs = fd.common_field((self.costs,costs),(1,1))
36		return Diagonal(self_costs*costs)
37
38	def norm(self,v):
39		costs,v = fd.common_field((self.costs,v),depths=(1,1))
40		return ad.Optimization.norm(costs*v,ord=2,axis=0)
41
42	def is_definite(self): return np.all(self.costs>0.,axis=0)
43	def anisotropy(self): return np.max(self.costs,axis=0)/np.min(self.costs,axis=0)
44	def cost_bound(self): return np.max(self.costs,axis=0)
45
46	@property
47	def vdim(self): return len(self.costs)
48	@property
49	def shape(self): return self.costs.shape[1:]
50
51	def flatten(self):      return self.costs
52	@classmethod
53	def expand(cls,arr):    return cls(arr)
54
55	def model_HFM(self):
56		return "Diagonal"+str(self.vdim)
57
58	@classmethod
59	def from_cast(cls,metric):
60		if isinstance(metric,cls): return metric
61		iso = Isotropic.from_cast(metric)
62		shape = (iso.vdim,) + iso.shape
63		return cls(np.broadcast_to(iso.cost,shape))
64
65	def __iter__(self):
66		yield self.costs
67
68	def make_proj_dual(self,newton_maxiter=7):
69		"""
70		Returns the orthogonal projection onto the unit ball.
71		- newton_maxiter : number of iterations for the inner subproblem
72		"""
73		# Implementation based on Benamou, J.-D., Carlier, G. & Hatchi, R. A numerical solution to 
74		# Monge’s problem with a Finsler distance as cost. ESAIM: M2AN 52, 2133–2148 (2018).
75		d_ = self.costs
76		def proj(x):
77			d,x = fd.common_field((d_,x),depths=(1,1))
78			x2=x**2
79			norm2 = np.sum(x2/d,axis=0)
80			# Solve, via a Newton method, the equation f(β)=0 where
81			# f(β) = -1 + sum_i xi^2/(1+di*β)^2 .
82			# By convexity of f, no damping or fancy stopping criterion is needed
83			β = np.min(d,axis=0)*(np.sqrt(norm2)-1) # Initial guess, exact in isotropic case.
84			orig = β<=0.; β[orig]=0. # Those are projected onto the origin
85			lc2 = np.where(orig[None],1.,d*x2)
86			for i in range(newton_maxiter):
87				ilb = 1./(d+β)
88				a = lc2*ilb**2
89				val = np.sum(a,axis=0)-1.
90				val[orig]=0
91				der = -2*np.sum(a*ilb,axis=0)
92				β -= val/der
93			return x*d/(d+β)
94		return proj

A Diagonal norm takes the form $$ F(x) = \sqrt{\sum_{0\leq i < d} c_i x_i^2 }, $$ where $(c_i)_{0\leq i < d}$, are given positive scalars

Member fields and __init__ arguments :

costs : the array of $(c_i)_{0\leq i < d}$ positive numbers. Required shape : $(d,n_1,..,n_k)$ where $d=$vdim is the ambient space dimension.

Diagonal(costs) View Source

23	def __init__(self,costs):
24		self.costs = ad.asarray(costs)

costs

@classmethod

def from_speed(cls, speeds): View Source

26	@classmethod
27	def from_speed(cls,speeds): 
28		"""
29		Produces a metric whose costs equal 1/speeds
30		"""
31		return cls(1./speeds)

Produces a metric whose costs equal 1/speeds

def dual(self): View Source

33	def dual(self): return self.from_speed(self.costs)

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

def with_costs(self, costs): View Source

34	def with_costs(self,costs): 
35		self_costs,costs = fd.common_field((self.costs,costs),(1,1))
36		return Diagonal(self_costs*costs)

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

def norm(self, v): View Source

38	def norm(self,v):
39		costs,v = fd.common_field((self.costs,v),depths=(1,1))
40		return ad.Optimization.norm(costs*v,ord=2,axis=0)

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 is_definite(self): View Source

42	def is_definite(self): return np.all(self.costs>0.,axis=0)

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

def anisotropy(self): View Source

43	def anisotropy(self): return np.max(self.costs,axis=0)/np.min(self.costs,axis=0)

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

def cost_bound(self): View Source

44	def cost_bound(self): return np.max(self.costs,axis=0)

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

vdim View Source

46	@property
47	def vdim(self): return len(self.costs)

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

shape View Source

48	@property
49	def shape(self): return self.costs.shape[1:]

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

def flatten(self): View Source

51	def flatten(self):      return self.costs

Flattens and concatenate the member fields into a single array.

@classmethod

def expand(cls, arr): View Source

52	@classmethod
53	def expand(cls,arr):    return cls(arr)

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

def model_HFM(self): View Source

55	def model_HFM(self):
56		return "Diagonal"+str(self.vdim)

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

@classmethod

def from_cast(cls, metric): View Source

58	@classmethod
59	def from_cast(cls,metric):
60		if isinstance(metric,cls): return metric
61		iso = Isotropic.from_cast(metric)
62		shape = (iso.vdim,) + iso.shape
63		return cls(np.broadcast_to(iso.cost,shape))

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

def make_proj_dual(self, newton_maxiter=7): View Source

68	def make_proj_dual(self,newton_maxiter=7):
69		"""
70		Returns the orthogonal projection onto the unit ball.
71		- newton_maxiter : number of iterations for the inner subproblem
72		"""
73		# Implementation based on Benamou, J.-D., Carlier, G. & Hatchi, R. A numerical solution to 
74		# Monge’s problem with a Finsler distance as cost. ESAIM: M2AN 52, 2133–2148 (2018).
75		d_ = self.costs
76		def proj(x):
77			d,x = fd.common_field((d_,x),depths=(1,1))
78			x2=x**2
79			norm2 = np.sum(x2/d,axis=0)
80			# Solve, via a Newton method, the equation f(β)=0 where
81			# f(β) = -1 + sum_i xi^2/(1+di*β)^2 .
82			# By convexity of f, no damping or fancy stopping criterion is needed
83			β = np.min(d,axis=0)*(np.sqrt(norm2)-1) # Initial guess, exact in isotropic case.
84			orig = β<=0.; β[orig]=0. # Those are projected onto the origin
85			lc2 = np.where(orig[None],1.,d*x2)
86			for i in range(newton_maxiter):
87				ilb = 1./(d+β)
88				a = lc2*ilb**2
89				val = np.sum(a,axis=0)-1.
90				val[orig]=0
91				der = -2*np.sum(a*ilb,axis=0)
92				β -= val/der
93			return x*d/(d+β)
94		return proj

Returns the orthogonal projection onto the unit ball.

newton_maxiter : number of iterations for the inner subproblem

Inherited Members

agd.Metrics.base.Base: gradient; disassociate; anisotropy_bound; cos_asym; cos; angle; inv_transform; 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