reg4opt package

Submodules

reg4opt.regression module

Solvers for operator and convex regression.

reg4opt.regression.closed_form_solution_cr(x, y, w, mu, L)

Closed form solution of convex regression with 2 training points.

Parameters
  • x (list) – The points \(x_i\) where the original function is evaluated.

  • y (list) – The evaluations of the function \(y_i = f(x_i)\).

  • w (list) – The valuations of the gradient \(w_i = \nabla f(x_i)\).

  • mu (float) – The strong convexity modulus of the target function.

  • L (float) – The smoothness modulus of the target function.

Returns

  • list – The approximate function evaluated in the points \(x_i\).

  • list – The gradient of the approximate function evaluated in the points \(x_i\).

reg4opt.regression.closed_form_solution_or(x, y, zeta)

Closed form solution of operator regression with 2 training points.

Parameters
  • x (list) – The points \(x_i\) where the original operator is evaluated.

  • y (list) – The evaluations of the operator \(y_i = \mathcal{T} x_i\).

  • zeta (float) – The Lipschitz modulus for the reconstructed operator.

Returns

The approximate operator evaluated in the points \(x_i\).

Return type

list

reg4opt.regression.convex_regression(x, y, mu, L, w=None, **kwargs)

Convex regression.

The function solves a convex regression problem

\[\begin{split}\begin{align} &\min_{f_i, g_i, \ i \in [D]} \frac{1}{D} \sum_{i = 1}^D (f_i - y_i)^2 \\ &\text{s.t.} \ f_i - f_j - \langle g_j, x_i - x_j \rangle \geq \\ &\qquad \frac{1}{2(1 - \mu / L)} \left( \frac{1}{L} \| g_i - g_j \|^2 + \mu \| x_i - x_j \|^2 - \frac{2\mu}{L} \langle g_j - g_i, x_j - x_i \rangle \right), \quad 1 \leq i < j \leq D \end{align}\end{split}\]

where \(y_i = f(x_i)\), \(i = 1, \ldots, D\), are evaluations of the function to be approximated, and \(0 < \mu < L\) are the strong convexity and smoothness moduli for the target function.

The function can also incorporate gradient information, in which case the problem becomes:

\[\begin{split}\begin{align} &\min_{f_i, g_i, \ i \in [D]} \frac{1}{D} \sum_{i = 1}^D (f_i - y_i)^2 + \| g_i - w_i \|^2 \\ &\text{s.t.} \ f_i - f_j - \langle g_j, x_i - x_j \rangle \geq \\ &\qquad \frac{1}{2(1 - \mu / L)} \left( \frac{1}{L} \| g_i - g_j \|^2 + \mu \| x_i - x_j \|^2 - \frac{2\mu}{L} \langle g_j - g_i, x_j - x_i \rangle \right), \quad 1 \leq i < j \leq D \end{align}\end{split}\]

with \(w_i = \nabla f_i(x_i)\).

In case \(D = 2\) (and gradient information is given), the function returns a closed form solution, see closed_form_solution_cr. Otherwise the solution is (approximately) computed using a solver based on the Peaceman-Rachford splitting, see prs_solver_cr.

Notice that the training data x and y (and the optional w) should be lists of column numpy arrays (for the sake of efficiency, no check is applied to them).

Parameters
  • x (list) – The points \(x_i\) where the original function is evaluated.

  • y (list) – The evaluations of the function \(y_i = f(x_i)\).

  • mu (float) – The strong convexity modulus of the target function.

  • L (float) – The smoothness modulus of the target function.

  • w (list, optional) – Optional gradient evaluations to add data to the regression problem.

  • **kwargs (tuple) – Arguments that should be passed to the solver.

Returns

  • list – The approximate function evaluated in the points \(x_i\).

  • list – The gradient of the approximate function evaluated in the points \(x_i\).

reg4opt.regression.operator_regression(x, y, zeta, **kwargs)

Operator regression.

The function solves an operator regression problem

\[\begin{split}\begin{align} &\min_{t_i, \ i \in [D]} \frac{1}{D} \sum_{i = 1}^D \| t_i - y_i \|^2 \\ &\text{s.t.} \ \| t_i - t_j \|^2 \leq \zeta^2 \| y_i - y_j \|^2, \quad 1 \leq i < j \leq D \end{align}\end{split}\]

where \(y_i = \mathcal{T} x_i\), \(i = 1, \ldots, D\), are evaluations of the operator to be approximated, and \(\zeta \in (0, 1]\) is the Lipschitz continuity modulus for the approximated operator.

In case \(D = 2\), the function returns a closed form solution, see closed_form_solution_or. If \(D > 2\) the solution is (approximately) computed using a solver based on the Peaceman-Rachford splitting, see prs_solver_or.

Notice that the training data x and y should be lists of column numpy arrays (for the sake of efficiency, no check is applied to them).

Parameters
  • x (list) – The points \(x_i\) where the original operator is evaluated.

  • y (list) – The evaluations of the operator \(y_i = \mathcal{T} x_i\).

  • zeta (float) – The Lipschitz modulus for the reconstructed operator.

  • **kwargs (tuple) – Optional arguments for the PRS solver.

Returns

  • list – The approximate operator evaluated in the points \(x_i\).

  • float – The Lipschitz modulus (either the given one or the one computed via autotuning).

reg4opt.regression.prs_solver_cr(x, y, mu, L, w=None, tol=0.01, num_iter=1000.0, rho=1)

Convex regression solver using PRS.

This function implements a tailored solver for convex regression based on the Peaceman-Rachford splitting. The solver can optionally exploit gradient information in the cost.

Parameters
  • x (list) – The points \(x_i\) where the original function is evaluated.

  • y (list) – The evaluations of the function \(y_i = f(x_i)\).

  • mu (float) – The strong convexity modulus of the target function.

  • L (float) – The smoothness modulus of the target function.

  • w (list, optional) – Optional evaluations of the gradient \(w_i = \nabla f(x_i)\).

  • tol (float, optional) – The solver stops if the fixed point residual \(\| x^{\ell+1} - x^\ell \|\) is below this tolerance.

  • num_iter (int, optional) – The maximum number of iterations that the solver can perform.

  • rho (float, optional) – The penalty parameter for PRS, must be a positive scalar.

Returns

  • list – The approximate function evaluated in the points \(x_i\).

  • list – The gradient of the approximate function evaluated in the points \(x_i\).

reg4opt.regression.prs_solver_or(x, y, zeta, tol=0.01, num_iter=1000.0, rho=1)

Operator regression solver using PRS.

This function implements a tailored solver for operator regression based on the Peaceman-Rachford splitting.

Parameters
  • x (list) – The points \(x_i\) where the original operator is evaluated.

  • y (list) – The evaluations of the operator \(y_i = \mathcal{T} x_i\).

  • zeta (float) – The Lipschitz modulus for the reconstructed operator.

  • tol (float, optional) – The solver stops if the fixed point residual \(\| x^{\ell+1} - x^\ell \|\) is below this tolerance.

  • num_iter (int, optional) – The maximum number of iterations that the solver can perform.

  • rho (float, optional) – The penalty parameter for PRS, must be a positive scalar.

Returns

The approximate operator evaluated in the points \(x_i\).

Return type

list

reg4opt.operators module

Operator template and examples.

class reg4opt.operators.AveragedOperator(o, a)

Bases: reg4opt.operators.Operator

Averaged operator.

This class defines an operator as the averaging (or relaxation) of a given operator. That is, given the operator \(\mathcal{T}\), and \(a \in (0,1]\), the class defines \((1 - a) \mathcal{I} + a \mathcal{T}\).

operator(x, *args, **kwargs)

An evaluation of the operator.

Parameters
  • x (array_like) – The x where the operator should be evaluated.

  • *args – The time at which the operator should be evaluated. Not required if the operator is static.

  • **kwargs – Any other required argument.

class reg4opt.operators.CompositeOperator(o2, o1)

Bases: reg4opt.operators.Operator

Composition of operators.

This class defines an operator as the composition of two given operators. That is, given the operators \(\mathcal{T}\) and \(\mathcal{R}\), the class defines \(\mathcal{T} \circ \mathcal{R}\). The domains and ranges of the two operators must be compatible.

operator(x, *args, **kwargs)

An evaluation of the operator.

Parameters
  • x (array_like) – The x where the operator should be evaluated.

  • *args – The time at which the operator should be evaluated. Not required if the operator is static.

  • **kwargs – Any other required argument.

class reg4opt.operators.DiscreteDynamicOperator(ops, t_s=1)

Bases: reg4opt.operators.Operator

Dynamic operator from a sequence of static ones.

This class creates a dynamic operator from a list of static operators. That is, given a sampling time \(T_\mathrm{s}\), the operator at time \(t_k = k T_\mathrm{s}\) is:

\[\mathcal{T}(\pmb{x}; t_k) = \mathcal{T}_k(\pmb{x})\]

with \(\mathcal{T}_k\) the k-th static operator in the list.

operator(x, t, **kwargs)

An evaluation of the operator.

Parameters
  • x (array_like) – The x where the operator should be evaluated.

  • *args – The time at which the operator should be evaluated. Not required if the operator is static.

  • **kwargs – Any other required argument.

sample(t)

Sample the operator.

The difference with the default Operator method is that it returns an operator in the list rather than a SampledOperator.

Parameters

t (float) – The time at which the operator should be sampled.

Returns

The closest operator in the list.

Return type

Operator

class reg4opt.operators.ForwardBackward(f, g, p)

Bases: reg4opt.operators.Operator

Forward-backward operator.

This class defines the forward-backward operator of two given cost functions \(f\) and \(g\). That is, the operator:

\[\operatorname{prox}_{\rho g}\left( x - \rho \nabla f(x) \right).\]
class reg4opt.operators.Gradient(f, s)

Bases: reg4opt.operators.Operator

Gradient operator.

This class defines the gradient step operator (or forward operator) of a given cost function \(f\).

operator(x, *args)

An evaluation of the operator.

Parameters
  • x (array_like) – The x where the operator should be evaluated.

  • *args – The time at which the operator should be evaluated. Not required if the operator is static.

  • **kwargs – Any other required argument.

class reg4opt.operators.Identity(dom, rng=None, time=None)

Bases: reg4opt.operators.Operator

Identity operator.

operator(x, *args, **kwargs)

An evaluation of the operator.

Parameters
  • x (array_like) – The x where the operator should be evaluated.

  • *args – The time at which the operator should be evaluated. Not required if the operator is static.

  • **kwargs – Any other required argument.

class reg4opt.operators.Operator(dom, rng=None, time=None)

Bases: object

Template class for operators.

This class serves as a template for operators:

\[\mathcal{T} : \mathbb{D}_1 \to \mathbb{D}_2\]

where \(\mathbb{D}_i \subset \mathbb{R}^{n_{i,1} \times n_{i,2} \times \cdots}\) for some given dimensions \(n_{i,1}, n_{i,2}, \ldots\), \(i = 1, 2\).

Operator objects support the following operations:

  • composition,

  • averaging (a.k.a. relaxation),

  • addition,

  • product by a scalar.

An Operator should expose the operator method to perform evaluations of \(\mathcal{T}\). The class then implements the fpr which computes the fixed point residual \(\| \mathcal{T} x - x\|\) at a given point \(x\).

The Operator can also be time-varying, in which case the function sample is also available.

dom

The operator’s domain.

Type

tvopt.sets.Set

rng

The operator’s range.

Type

tvopt.sets.Set

time

The time domain \(\mathbb{R}_+\). If the operator is static this is None.

Type

tvopt.sets.T

is_dynamic

Attribute to check if the operator is static or dynamic.

Type

bool

average(a)

Averaging (a.k.a. relaxation) of the operator.

The method returns an averaged version of the operator, defined as

\[(1 - a) \mathcal{I} + a \mathcal{T}\]

where \(a \in (0,1]\) and \(\mathcal{I}\) is the identity operator.

Parameters

a (float) – The relaxation constant.

Returns

DESCRIPTION.

Return type

Operator

See also

AveragedOperator

The averaged operator object.

compose(other)

Compose the operator with a second one.

The method composes the operator with a second one, provided that the range of the former coincides with the domain of the latter.

Parameters

other (Operator) – The second operator.

Returns

An operator object defining the composition.

Return type

Operator

See also

CompositeOperator

The composite operator object.

fpr(x, *args, **kwargs)

Evaluate the fixed point residual (FPR).

This method returns the FPR

\[\| \mathcal{T} x - x \|.\]
Parameters
  • x (array_like) – The point whose FPR should be evaluated.

  • *args – The time at which the operator should be evaluated. Not required if the operator is static.

  • **kwargs – Any other required argument.

Returns

The FPR at x.

Return type

float

operator(x, *args, **kwargs)

An evaluation of the operator.

Parameters
  • x (array_like) – The x where the operator should be evaluated.

  • *args – The time at which the operator should be evaluated. Not required if the operator is static.

  • **kwargs – Any other required argument.

sample(t)

Sample the operator.

This method returns a SampledOperator object which consists of the operator as evaluated at time t.

If the operator is static, the operator itself is returned.

Parameters

t (float) – The time at which the operator should be sampled.

Returns

The sampled operator or, if static, the operator itself.

Return type

Operator

See also

SampledOperator

The sampled operator object.

class reg4opt.operators.PeacemanRachford(f, g, p)

Bases: reg4opt.operators.Operator

Peaceman-Rachford operator.

This class defines the Peaceman-Rachford operator of two given cost functions \(f\) and \(g\). That is, the operator:

\[\operatorname{refl}_{\rho g} \left( \operatorname{refl}_{\rho f}(x) \right).\]
class reg4opt.operators.Proximal(f, p)

Bases: reg4opt.operators.Operator

Proximal operator.

This class defines the proximal operator (or backward operator) of a given cost function \(f\). That is, the operator:

\[\operatorname{prox}_{\rho f}(x) = \operatorname{arg\,min}_{y} \left\{ f(y) + \frac{1}{2 \rho} \|y - x\|^2 \right\}.\]
operator(x, *args)

An evaluation of the operator.

Parameters
  • x (array_like) – The x where the operator should be evaluated.

  • *args – The time at which the operator should be evaluated. Not required if the operator is static.

  • **kwargs – Any other required argument.

class reg4opt.operators.Reflective(f, p)

Bases: reg4opt.operators.Operator

Reflective operator.

This class defines the reflective operator of a given cost function \(f\). That is, the operator:

\[\operatorname{refl}_{\rho f}(x) = 2 \operatorname{prox}_{\rho f}(x) - x.\]
class reg4opt.operators.SampledOperator(o, t)

Bases: reg4opt.operators.Operator

Sampled operator.

This class defines a static operator by sampling a given dynamic operator, that is, by fixing the time argument to a given value.

operator(x, **kwargs)

An evaluation of the operator.

Parameters
  • x (array_like) – The x where the operator should be evaluated.

  • *args – The time at which the operator should be evaluated. Not required if the operator is static.

  • **kwargs – Any other required argument.

class reg4opt.operators.ScaledOperator(o, c)

Bases: reg4opt.operators.Operator

Scaled operator.

This class defines an operator multiplied by a scalar. That is, given the operator \(\mathcal{T}\) and \(c \in \mathbb{R}\), the class defines \(c \mathcal{T}\).

operator(x, *args, **kwargs)

An evaluation of the operator.

Parameters
  • x (array_like) – The x where the operator should be evaluated.

  • *args – The time at which the operator should be evaluated. Not required if the operator is static.

  • **kwargs – Any other required argument.

class reg4opt.operators.SeparableOperator(ops)

Bases: reg4opt.operators.Operator

Separable operator.

This class defines a separable operator, that is

\[\begin{split}\mathcal{T}(\pmb{x}; t) = \begin{bmatrix} \vdots \\ \mathcal{T}_i(x_i; t) \\ \vdots \end{bmatrix}\end{split}\]

where \(x_i \in \mathbb{R}^{n_1 \times n_2 \times \ldots}\) for each \(i = 1, \ldots, N\). Each of the component operators \(T_i\) can be either static or dynamic. This is useful for defining distributed optimization problems.

The overall dimension of the domain is \(n_1 \times n_2 \times \ldots \times N\), meaning that the last dimension indexes the components.

The class exposes the same methods as any Operator, with the difference that the keyword argument i can be used to evaluate only a single component. If all components are evaluated, an ndarray is returned with the last dimension indexing the components.

The class has the Operator attributes, with the following additions.

ops

The component operators.

Type

list

N

The number of components.

Type

int

is_dynamic

True if at least one component is dynamic.

Type

bool

operator(x, *args, i=None, **kwargs)

An evaluation of the separable operator.

This method performs an evaluation of the component operators. If the keyword argument i is specified, then only the corresponding operator is evaluated.

Parameters
  • x (array_like) – The x where the operator(s) should be evaluated.

  • *args – The time at which the operator(s) should be evaluated. Not required if they are static.

  • i (int, optional) – If specified, only the corresponding component operator \(\mathcal{T}_i\) is evaluated.

  • **kwargs – Any other required argument.

Returns

If i is specified, the evaluation of the i-th component, otherwise an ndarray stacking the components evaluations along the last dimension.

Return type

ndarray

class reg4opt.operators.SumOperator(o1, o2)

Bases: reg4opt.operators.Operator

Sum of operators.

This class defines an operator as the sum of two given operators. That is, given the operators \(\mathcal{T}\) and \(\mathcal{R}\), the class defines \(\mathcal{T} + \mathcal{R}\). The domains and ranges of the two operators must be compatible.

operator(x, *args, **kwargs)

An evaluation of the operator.

Parameters
  • x (array_like) – The x where the operator should be evaluated.

  • *args – The time at which the operator should be evaluated. Not required if the operator is static.

  • **kwargs – Any other required argument.

reg4opt.interpolation module

Interpolation of Lipschitz continuous operators.

reg4opt.interpolation.alternating_projections(s, x=None, tol=1e-10, num_iters=1000)

Method of alternating projections.

Given an initial point \(x^0\), this method computes for \(\ell \in \mathbb{N}\):

\[x^{\ell+1} = \operatorname{proj}_n \circ \cdots \circ \operatorname{proj}_1 x^\ell\]

where \(\operatorname{proj}_i\) is the projection operator onto the \(i\)-th set in the list s. See 1 for details and convergence.

Parameters
  • s (list) – A list of Set objects defining the problem.

  • x (array_like, optional) – The initial point \(x^0\), defaults to a randomly picked one if not given.

  • tol (float, optional) – The algorithm is stopped if the fixed point residual \(\|x^{\ell+1} - x^\ell\|\) is smaller than tol.

  • num_iters (int, optional) – The maximum number of iterations that the algorithm can perform.

Returns

x – A point in the intersection of the given sets.

Return type

ndarray

References

1

H. Bauschke and V. Koch, “Projection Methods: Swiss Army Knives for Solving Feasibility and Best Approximation Problems with Halfspaces,” in Contemporary Mathematics, vol. 636, S. Reich and A. Zaslavski, Eds. Providence, Rhode Island: American Mathematical Society, 2015, pp. 1–40.

reg4opt.interpolation.dykstra(s, v, tol=1e-10, num_iters=1000)

Dykstra projection method.

Given an initial point \(x^0 = v\), this method computes for \(\ell \in \mathbb{N}\):

\[\begin{split}\begin{align} \bar{x}^{\ell+1} &= \frac{1}{n} \sum_{i = 1}^n (x_i^\ell + p_i^\ell) \\ p_i^{\ell+1} &= x_i^\ell + p_i^\ell - \bar{x}^{\ell+1} \\ x_i^{\ell+1} &= \operatorname{proj}_i(\bar{x}^{\ell+1} + q_i^\ell) \\ q_i^{\ell+1} &= \bar{x}^{\ell+1} + q_i^\ell - x_i^{\ell+1} \end{align}\end{split}\]

where \(\operatorname{proj}_i\) is the projection operator onto the \(i\)-th set in the list s. See 2 and 3 for details and convergence.

Parameters
  • s (list) – A list of Set objects defining the problem.

  • v (array_like) – The point to be projected.

  • tol (float, optional) – The algorithm is stopped if the fixed point residual \(\|x^{\ell+1} - x^\ell\|\) is smaller than tol.

  • num_iters (int, optional) – The maximum number of iterations that the algorithm can perform.

Returns

x – An approximate projection of v onto the intersection of the sets.

Return type

ndarray

References

2

H. Bauschke and V. Koch, “Projection Methods: Swiss Army Knives for Solving Feasibility and Best Approximation Problems with Halfspaces,” in Contemporary Mathematics, vol. 636, S. Reich and A. Zaslavski, Eds. Providence, Rhode Island: American Mathematical Society, 2015, pp. 1–40.

3

H. H. Bauschke, R. S. Burachik, D. B. Herman, and C. Y. Kaya, “On Dykstra’s algorithm: finite convergence, stalling, and the method of alternating projections,” Optimization Letters, vol. 14, no. 8, pp. 1975–1987, Nov. 2020.

reg4opt.interpolation.halpern(s, v, tol=1e-10, num_iters=1000)

Halpern projection method.

Given an initial point \(x^0 = v\), this method computes for \(\ell \in \mathbb{N}\):

\[\begin{split}\begin{align} y^{\ell+1} &= \operatorname{proj}_n \circ \cdots \circ \operatorname{proj}_1 x^\ell \\ x^{\ell+1} &= \frac{v}{l+1} + \frac{l}{l+1} y^{\ell+1} \end{align}\end{split}\]

where \(\operatorname{proj}_i\) is the projection operator onto the \(i\)-th set in the list s. See 4 and 5 for details and convergence.

Parameters
  • s (list) – A list of Set objects defining the problem.

  • v (array_like) – The point to be projected.

  • tol (float, optional) – The algorithm is stopped if the fixed point residual \(\|x^{\ell+1} - x^\ell\|\) is smaller than tol.

  • num_iters (int, optional) – The maximum number of iterations that the algorithm can perform.

Returns

x – An approximate projection of v onto the intersection of the sets.

Return type

ndarray

References

4

H. Bauschke and V. Koch, “Projection Methods: Swiss Army Knives for Solving Feasibility and Best Approximation Problems with Halfspaces,” in Contemporary Mathematics, vol. 636, S. Reich and A. Zaslavski, Eds. Providence, Rhode Island: American Mathematical Society, 2015, pp. 1–40.

5

H. H. Bauschke and P. L. Combettes, Convex analysis and monotone operator theory in Hilbert spaces, 2nd ed. Cham: Springer, 2017.

reg4opt.interpolation.haugazeau(s, v, tol=1e-10, num_iters=1000)

Haugazeau projection method.

Given an initial point \(x^0 = v\), this method computes for \(\ell \in \mathbb{N}\):

\[x^{\ell+1} = Q(v, x^\ell, \operatorname{proj}_i x^\ell)\]

choosing \(i = 1, \ldots, n\) sequentially, and where \(\operatorname{proj}_i\) is the projection operator onto the \(i\)-th set in the list s, and

\[\begin{split}Q(x, y, z) = \begin{cases} z & \text{if} \ d = 0, a \geq 0 \\ x + (1 + a/c) (z - y) & \text{if} \ d > 0, a c \geq d \\ y + (c/d) (a(x - y) + b(z - y)) & \text{if} \ d > 0, a c < d \\ \text{undefined otherwise} \end{cases}\end{split}\]

where \(a = (x - y)^\top (y - z)\), \(b = \|x - y\|^2\), \(c = \|y - z\|^2\), and \(d = b c - a^2\). See 6 for details and convergence.

Parameters
  • s (list) – A list of Set objects defining the problem.

  • v (array_like) – The point to be projected.

  • tol (float, optional) – The algorithm is stopped if the fixed point residual \(\|x^{\ell+1} - x^\ell\|\) is smaller than tol.

  • num_iters (int, optional) – The maximum number of iterations that the algorithm can perform.

Returns

x – An approximate projection of v onto the intersection of the sets.

Return type

ndarray

References

6

H. Bauschke and V. Koch, “Projection Methods: Swiss Army Knives for Solving Feasibility and Best Approximation Problems with Halfspaces,” in Contemporary Mathematics, vol. 636, S. Reich and A. Zaslavski, Eds. Providence, Rhode Island: American Mathematical Society, 2015, pp. 1–40.

reg4opt.interpolation.interpolator(x, x_i, t_i, zeta, t0=None, solver=None, **solver_params)

Operator interpolation.

This function interpolates the \(\zeta\)-Lipschitz continuous operator \(\mathcal{T} : \mathbb{R}^n \to \mathbb{R}^n\) to the point \(x\), preserving the Lipschitz continuity.

In particular, given the operator evaluations

\[t_i = \mathcal{T} x_i, \quad i = 1, \ldots, D\]

the function computes \(\hat{\mathcal{T}} x\) as

\[\begin{split}\hat{\mathcal{T}} x = \begin{cases} t_i & \text{if} \ x = x_i \\ \hat{t} \in \bigcap_{i \in [D]} \mathbb{B}_{\zeta \|x - x_i\|}(t_i) & \text{otherwise} \end{cases}\end{split}\]

where \(\mathbb{B}_{\zeta \|x - x_i\|}(t_i)\) is a ball of center \(t_i\) and radius \(\zeta \|x - x_i\|\). We see then that interpolating requires finding a point in the intersection of balls. To this end, an efficient approach is to use the method of alternating projections (MAP) (see alternating_projections), but other solvers are available.

Parameters
  • x (array_like) – The point where the interpolated operator should be evaluated.

  • x_i (list) – The points of the regression problem.

  • t_i (list) – The solution of the regression problem.

  • zeta (float) – The contraction constant.

  • t0 (array_like, optional) – The initial condition for the projection algorithm, defaults to a vector of zeros.

  • solver (str) – The method that should be used to find a point in the intersection of balls. By default the MAP is used, due to its efficiency in this scenario.

  • solver_params (tuple) – Parameters for the solver, for example the maximum number of iterations that it can perform.

Returns

The interpolated operator evaluated at x.

Return type

ndarray

See also

alternating_projections

The method of alternating projections (MAP).

Notes

If the operator is defined in \(\mathbb{R}\) (that is, \(n = 1\)), then the balls become closed intervals on the line and the intersection can be computed in closed form.

reg4opt.interpolation.parallel_projections(s, x=None, tol=1e-10, num_iters=1000)

Method of parallel projections.

Given an initial point \(x^0\), this method computes for \(\ell \in \mathbb{N}\):

\[x^{\ell+1} = \frac{1}{n} \sum_{i = 1}^n \operatorname{proj}_i x^\ell\]

where \(\operatorname{proj}_i\) is the projection operator onto the \(i\)-th set in the list s. See 7 for details and convergence.

Parameters
  • s (list) – A list of Set objects defining the problem.

  • x (array_like, optional) – The initial point \(x^0\), defaults to a randomly picked one if not given.

  • tol (float, optional) – The algorithm is stopped if the fixed point residual \(\|x^{\ell+1} - x^\ell\|\) is smaller than tol.

  • num_iters (int, optional) – The maximum number of iterations that the algorithm can perform.

Returns

x – A point in the intersection of the given sets.

Return type

ndarray

References

7

H. Bauschke and V. Koch, “Projection Methods: Swiss Army Knives for Solving Feasibility and Best Approximation Problems with Halfspaces,” in Contemporary Mathematics, vol. 636, S. Reich and A. Zaslavski, Eds. Providence, Rhode Island: American Mathematical Society, 2015, pp. 1–40.

reg4opt.interpolation.peaceman_rachford(s, x=None, a=0.95, tol=1e-10, num_iters=1000)

Peaceman-Rachford splitting (PRS).

Given an initial point \(x^0\), this method computes for \(\ell \in \mathbb{N}\):

\[\begin{split}\begin{align} x_i^{\ell+1} &= \operatorname{proj}_i z_i^\ell \\ \bar{x}^{\ell+1} &= \frac{1}{n} \sum_{i = 1}^n x_i^{\ell+1} \\ z_i^{\ell+1} &= (1-a) z_i^\ell + a (2 \bar{x}^{\ell+1} - x_i^{\ell+1}) \end{align}\end{split}\]

where \(\operatorname{proj}_i\) is the projection operator onto the \(i\)-th set in the list s, and \(a \in (0,1]\). See 8 for details and convergence.

Parameters
  • s (list) – A list of Set objects defining the problem.

  • x (array_like, optional) – The initial point \(x^0\), defaults to a randomly picked one if not given.

  • a (float, optional) – The relaxation constant of the Peaceman-Rachford, defaults to 0.95.

  • tol (float, optional) – The algorithm is stopped if the fixed point residual \(\|x^{\ell+1} - x^\ell\|\) is smaller than tol.

  • num_iters (int, optional) – The maximum number of iterations that the algorithm can perform.

Returns

x – A point in the intersection of the given sets.

Return type

ndarray

References

8

H. Bauschke and V. Koch, “Projection Methods: Swiss Army Knives for Solving Feasibility and Best Approximation Problems with Halfspaces,” in Contemporary Mathematics, vol. 636, S. Reich and A. Zaslavski, Eds. Providence, Rhode Island: American Mathematical Society, 2015, pp. 1–40.

reg4opt.utils module

Utility functions.

reg4opt.utils.generate_data(T, x, num_data, method='normal', **kwargs)

Generate the training data for an operator regression problem.

The function encapsulates different method of generating the training data for an operator regression problem.

Parameters
  • T (operators.Operator) – The operator source of the training data.

  • x (ndarray) – The center of the training data, that is, all other training points are chosen as perturbations of x.

  • num_data (int) – The number of training points, including x, which means that num_data-1 points are randomly selected.

  • method (str, optional) – The method to choose the training points, defaults to “normal”.

  • **kwargs (tuple) – The arguments of method.

Returns

  • list – The points \(x_i\) selected by the method.

  • list – The operator T evaluated in the chosen points.

See also

sample_training_points

The function used to sample \(x_i\).

reg4opt.utils.generate_data_cr(f, x, num_data, gradient=False, method='normal', **kwargs)

Generate the training data for a convex regression problem.

The function encapsulates different method of generating the training data for a convex regression problem, including gradient information if required.

Parameters
  • f (tvopt.costs.Cost) – The cost source of the training data.

  • x (ndarray) – The center of the training data, that is, all other training points are chosen as perturbations of x.

  • num_data (int) – The number of training points, including x, which means that num_data-1 points are randomly selected.

  • gradient (bool, optional) – Specify if also gradient information should be returned.

  • method (str, optional) – The method to choose the training points, defaults to “normal”.

  • **kwargs (tuple) – The arguments of method.

Returns

  • list – The points \(x_i\) selected by the method.

  • list – The operator T evaluated in the chosen points.

See also

sample_training_points

The function used to sample \(x_i\).

reg4opt.utils.norm(x)

Compute the norm of the given vector.

Parameters

x (array_like) – The vector array.

Returns

The square norm.

Return type

ndarray

See also

square_norm

Square norm

Notes

The function reshapes x to a column vector, so it does not correctly handle n-dimensional arrays. For n-dim arrays use numpy.linalg.norm.

reg4opt.utils.norm_1(x)

Compute the \(\ell_1\) norm of the given vector.

Parameters

x (array_like) – The vector array.

Returns

The \(\ell_1\) norm.

Return type

ndarray

See also

norm

Euclidean norm

Notes

The function treats x as a one-dimensional vector (row or column). For n-dim arrays use numpy.linalg.norm.

reg4opt.utils.print_progress(i, num_iter, bar_length=80, decimals=2)

Print the progresso to command line.

Parameters
  • i (int) – Current iteration.

  • num_iter (int) – Total number of iterations.

  • bar_length (int, optional) – Length of progress bar.

  • decimals (int, optional) – Decimal places of the progress percent.

Notes

Adapted from here.

reg4opt.utils.random_intersecting_balls(n, d, max_radius=100)

Randomly generate a set of intersecting ball sets.

This utility function generates a set of d balls in \(\mathbb{R}^n\) for the purpose of testing projection methods. The centers of the balls are uniformly chosen points in the unit sphere, and the radius is randomly picked between 0 (excluded) and max_radius.

As a consequence, the origin of \(\mathbb{R}^n\) is (one of the) points in the intersection of the balls.

Parameters
  • n (int) – The dimension of the space.

  • d (int) – The number of balls to be generated.

  • max_radius (float, optional) – The maximum radius of the generated balls. The default is 100.

Returns

s – A list of tvopt.sets.Ball objects representing the sets.

Return type

list

Notes

Adapted from here.

reg4opt.utils.sample_training_points(x, num_data, method='normal', **kwargs)

Sample training points.

The function encapsulates different methods of generating the training points for a regression problem. The choices are

  • normal: the data are chosen as \(x + d_i\), where \(d_i\) are random vectors with normal distribution (by default the variance is 0.01);

  • fireworks: the data are chosen as \(x + d_i\), where we have \(d_i = a e_j\), with \(a\) a r.v. with normal distribution (with default variance 0.01) and \(e_j\) a randomly selected vector of the standard basis;

  • uniform: the data are chosen as \(x + d_i\), where \(d_i\) are random vectors with uniform distribution in [-a, a] (by default a = 1).

Parameters
  • x (ndarray) – The center of the training data, that is, all other training points are chosen as perturbations of x.

  • num_data (int) – The number of training points, including x, which means that num_data-1 points are randomly selected.

  • method (str, optional) – The method to choose the training points, defaults to “normal”. Alternatives are “fireworks” (alias “fw”, “f”) and “uniform” (alias “u”).

  • **kwargs (tuple) – The arguments of method.

Returns

The points \(x_i\) selected by the method.

Return type

list

reg4opt.utils.square_norm(x)

Compute the square norm of the given vector.

Parameters

x (array_like) – The vector array.

Returns

The square norm.

Return type

ndarray

Notes

The function reshapes x to a column vector, so it does not correctly handle n-dimensional arrays. For n-dim arrays use numpy.linalg.norm.

Module contents