import numpy as np
from pylops.utils.backend import get_module, to_numpy
[docs]def gradtest_proximal(Op, n, x=None, dtype="float64",
delta=1e-6, rtol=1e-6, atol=1e-21,
complexflag=False, raiseerror=True,
verb=False, backend="numpy"):
r"""Gradient test for Proximal operator.
Compute the gradient of ``Op`` using both the provided method and a
numerical approximation with a perturbation ``delta`` applied to a
single, randomly selected parameter of the input vector.
Parameters
----------
Op : :obj:`pyproximal.Proximal`
Proximal operator to test.
n : :obj:`int`
Size of input vector
x : :obj:`numpy.ndarray`, optional
Input vector (if ``None``, randomly drawn from a
Normal distribution)
dtype : :obj:`str`, optional
Dtype of vector ``x`` to generate (only used when ``x=None``)
delta : :obj:`float`, optional
Perturbation
rtol : :obj:`float`, optional
Relative gradtest tolerance
atol : :obj:`float`, optional
Absolute gradtest tolerance
complexflag : :obj:`bool`, optional
Generate random vectors with real (``False``) or
complex (``True``) entries
raiseerror : :obj:`bool`, optional
Raise error or simply return ``False`` when dottest fails
verb : :obj:`bool`, optional
Verbosity
backend : :obj:`str`, optional
Backend used for dot test computations (``numpy`` or ``cupy``). This
parameter will be used to choose how to create the random vectors.
Returns
-------
passed : :obj:`bool`
Passed flag.
Raises
------
AssertionError
If grad-test is not verified within chosen tolerances.
Notes
-----
A gradient-test is mathematical tool used in the development of numerical
nonliner operators.
More specifically, a correct implementation of the gradient for
a nonlinear operator should verify the following *equality*
within a numerical tolerance:
.. math::
\frac{\partial Op(\mathbf{x})}{\partial \mathbf{x}} =
\frac{Op(\mathbf{x}+\delta \mathbf{x})-Op(\mathbf{x})}{\delta \mathbf{x}}
"""
ncp = get_module(backend)
# get random vectors for x and y
if x is None:
x = np.random.normal(0., 1., n).astype(dtype)
if complexflag:
x = x + 1j * np.random.normal(0., 1., n).astype(dtype)
# compute function
f = Op(x)
# compute gradient
g = Op.grad(x)
# choose location of perturbation, whether to act on x or y and on real or imag part
iqx = np.random.randint(0, n)
r_or_i = np.random.randint(0, 2)
if r_or_i == 0:
delta1 = delta
else:
delta1 = delta * 1j
# extract gradient value to test
x[iqx] = x[iqx] + delta1
grad = g[iqx]
# compute new function at perturbed location
fdelta = Op(x)
# evaluate if gradient test passed
grad_delta = (fdelta - f) / np.abs(delta)
grad_diff = grad_delta - (grad.real if r_or_i == 0 else grad.imag)
passed = np.isclose(grad_diff, 0, rtol, atol)
# verbosity or error raising
if (not passed and raiseerror) or verb:
passed_status = "passed" if passed else "failed"
msg = f"Grad test {passed_status}, Analytic={grad.real if r_or_i == 0 else grad.imag} - " \
f"Numeric={grad_delta}"
if not passed and raiseerror:
raise AssertionError(msg)
else:
print(msg)
return passed
[docs]def gradtest_bilinear(Op, delta=1e-6, rtol=1e-6, atol=1e-21,
complexflag=False, raiseerror=True,
verb=False, backend="numpy"):
r"""Gradient test for Bilinear operator.
Compute the gradient of ``Op`` using both the provided method and a
numerical approximation with a perturbation ``delta`` applied to a
single, randomly selected parameter of either the ``x`` or ``y``
vectors.
Parameters
----------
Op : :obj:`pyproximal.utils.BilinearOperator`
Bilinear operator to test.
delta : :obj:`float`, optional
Perturbation
rtol : :obj:`float`, optional
Relative gradtest tolerance
atol : :obj:`float`, optional
Absolute gradtest tolerance
complexflag : :obj:`bool`, optional
Generate random vectors with real (``False``) or
complex (``True``) entries
raiseerror : :obj:`bool`, optional
Raise error or simply return ``False`` when dottest fails
verb : :obj:`bool`, optional
Verbosity
backend : :obj:`str`, optional
Backend used for dot test computations (``numpy`` or ``cupy``). This
parameter will be used to choose how to create the random vectors.
Returns
-------
passed : :obj:`bool`
Passed flag.
Raises
------
AssertionError
If grad-test is not verified within chosen tolerances.
Notes
-----
A gradient-test is mathematical tool used in the development of numerical
bilinear operators.
More specifically, a correct implementation of the gradient for
a bilinear operator should verify the following *equalities*
within a numerical tolerance:
.. math::
\frac{\partial Op(\mathbf{x})}{\partial \mathbf{x}} =
\frac{Op(\mathbf{x}+\delta \mathbf{x}, \mathbf{y})-
Op(\mathbf{x})}{\delta \mathbf{x}, \mathbf{y}}
and
.. math::
\frac{\partial Op(\mathbf{x}, \mathbf{y})}{\partial \mathbf{y}} =
\frac{Op(\mathbf{x}, \mathbf{y}+\delta \mathbf{y})-
Op(\mathbf{x}, \mathbf{y})}{\delta \mathbf{y}}
"""
ncp = get_module(backend)
nx = Op.sizex
ny = Op.sizey
# extract x and y from Op
x, y = Op.x.ravel(), Op.y.ravel()
# compute function at x and y
f = Op(x, y)
# compute gradients at x and y
gx = Op.gradx(x)
gy = Op.grady(y)
# choose location of perturbation, whether to act on x or y and on real or imag part
iqx, iqy = np.random.randint(0, nx), np.random.randint(0, ny)
x_or_y = np.random.randint(0, 2)
delta1 = delta
if complexflag:
r_or_i = np.random.randint(0, 2)
if r_or_i == 1:
delta1 = delta * 1j
# extract gradient value to test
if x_or_y == 0:
x[iqx] = x[iqx] + delta1
grad = gx[iqx]
else:
y[iqy] = y[iqy] + delta1
grad = gy[iqy]
# compute new function at perturbed location
fdelta = Op(x, y)
# evaluate if gradient test passed
grad_delta = (fdelta - f) / np.abs(delta)
grad_diff = grad_delta - (grad.real if not complexflag or r_or_i == 0 else grad.imag)
passed = np.isclose(grad_diff, 0, rtol, atol)
# verbosity or error raising
if (not passed and raiseerror) or verb:
passed_status = "passed" if passed else "failed"
msg = f"Grad test {passed_status}, Analytic={grad.real if r_or_i == 0 else grad.imag} - " \
f"Numeric={grad_delta}"
if not passed and raiseerror:
raise AssertionError(msg)
else:
print(msg)
return passed