"""
Functions concerned with solving for decoders or full weight matrices.
Many of the solvers in this file can solve for decoders or weight matrices,
depending on whether the post-population encoders `E` are provided (see below).
Solvers that are only intended to solve for either decoders or weights can
remove the `E` parameter or make it manditory as they see fit.
"""
import collections
import logging
import time
import numpy as np
import nengo.utils.least_squares_solvers as lstsq
from nengo.exceptions import ValidationError
from nengo.params import BoolParam, NumberParam, Parameter
from nengo.utils.compat import range, with_metaclass, iteritems
from nengo.utils.least_squares_solvers import (
format_system, rmses, LeastSquaresSolverParam)
from nengo.utils.magic import DocstringInheritor
logger = logging.getLogger(__name__)
[docs]class Solver(with_metaclass(DocstringInheritor)):
"""Decoder or weight solver."""
[docs] def __call__(self, A, Y, rng=None, E=None):
"""Call the solver.
Parameters
----------
A : (n_eval_points, n_neurons) array_like
Matrix of the neurons' activities at the evaluation points
Y : (n_eval_points, dimensions) array_like
Matrix of the target decoded values for each of the D dimensions,
at each of the evaluation points.
rng : `numpy.random.RandomState`, optional (Default: None)
A random number generator to use as required. If None,
the ``numpy.random`` module functions will be used.
E : (dimensions, post.n_neurons) array_like, optional (Default: None)
Array of post-population encoders. Providing this tells the solver
to return an array of connection weights rather than decoders.
Returns
-------
X : (n_neurons, dimensions) or (n_neurons, post.n_neurons) ndarray
(n_neurons, dimensions) array of decoders (if ``solver.weights``
is False) or (n_neurons, post.n_neurons) array of weights
(if ``'solver.weights`` is True).
info : dict
A dictionary of information about the solver. All dictionaries have
an ``'rmses'`` key that contains RMS errors of the solve.
Other keys are unique to particular solvers.
"""
raise NotImplementedError("Solvers must implement '__call__'")
def __hash__(self):
items = list(self.__dict__.items())
items.sort(key=lambda item: item[0])
hashes = []
for k, v in items:
if isinstance(v, np.ndarray):
if v.size < 1e5:
a = v[:]
a.setflags(write=False)
hashes.append(hash(a))
else:
raise ValidationError("array is too large to hash", attr=k)
elif isinstance(v, collections.Iterable):
hashes.append(hash(tuple(v)))
elif isinstance(v, collections.Callable):
hashes.append(hash(v.__code__))
else:
hashes.append(hash(v))
return hash(tuple(hashes))
def __str__(self):
return "%s(%s)" % (
self.__class__.__name__,
', '.join("%s=%s" % (k, v) for k, v in iteritems(self.__dict__)))
[docs] def mul_encoders(self, Y, E, copy=False):
"""Helper function that projects signal ``Y`` onto encoders ``E``.
Parameters
----------
Y : ndarray
The signal of interest.
E : (dimensions, n_neurons) array_like or None
Array of encoders. If None, ``Y`` will be returned unchanged.
copy : bool, optional (Default: False)
Whether a copy of ``Y`` should be returned if ``E`` is None.
"""
if self.weights:
if E is None:
raise ValidationError(
"Encoders must be provided for weight solver", attr='E')
return np.dot(Y, E)
else:
if E is not None:
raise ValidationError(
"Encoders must be 'None' for decoder solver", attr='E')
return Y.copy() if copy else Y
class SolverParam(Parameter):
def validate(self, instance, solver):
if solver is not None and not isinstance(solver, Solver):
raise ValidationError("'%s' is not a solver" % solver,
attr=self.name, obj=instance)
super(SolverParam, self).validate(instance, solver)
[docs]class Lstsq(Solver):
"""Unregularized least-squares solver.
Parameters
----------
weights : bool, optional (Default: False)
If False, solve for decoders. If True, solve for weights.
rcond : float, optional (Default: 0.01)
Cut-off ratio for small singular values (see `numpy.linalg.lstsq`).
Attributes
----------
rcond : float
Cut-off ratio for small singular values (see `numpy.linalg.lstsq`).
weights : bool
If False, solve for decoders. If True, solve for weights.
"""
def __init__(self, weights=False, rcond=0.01):
self.rcond = rcond
self.weights = weights
def __call__(self, A, Y, rng=None, E=None):
tstart = time.time()
Y = self.mul_encoders(Y, E)
X, residuals2, rank, s = np.linalg.lstsq(A, Y, rcond=self.rcond)
t = time.time() - tstart
return X, {'rmses': rmses(A, X, Y),
'residuals': np.sqrt(residuals2),
'rank': rank,
'singular_values': s,
'time': t}
class _LstsqNoiseSolver(Solver):
"""Base class for least-squares solvers with noise."""
weights = BoolParam('weights')
noise = NumberParam('noise', low=0)
solver = LeastSquaresSolverParam('solver')
def __init__(self, weights=False, noise=0.1, solver=lstsq.Cholesky()):
"""
Parameters
----------
weights : bool, optional (Default: False)
If False, solve for decoders. If True, solve for weights.
noise : float, optional (Default: 0.1)
Amount of noise, as a fraction of the neuron activity.
solver : `.LeastSquaresSolver`, optional (Default: ``Cholesky()``)
Subsolver to use for solving the least squares problem.
Attributes
----------
noise : float
Amount of noise, as a fraction of the neuron activity.
solver : `.LeastSquaresSolver`
Subsolver to use for solving the least squares problem.
weights : bool
If False, solve for decoders. If True, solve for weights.
"""
self.weights = weights
self.noise = noise
self.solver = solver
[docs]class LstsqNoise(_LstsqNoiseSolver):
"""Least-squares solver with additive Gaussian white noise."""
def __call__(self, A, Y, rng=None, E=None):
tstart = time.time()
rng = np.random if rng is None else rng
sigma = self.noise * A.max()
A = A + rng.normal(scale=sigma, size=A.shape)
X, info = self.solver(A, Y, 0, rng=rng)
info['time'] = time.time() - tstart
return self.mul_encoders(X, E), info
[docs]class LstsqMultNoise(_LstsqNoiseSolver):
"""Least-squares solver with multiplicative white noise."""
def __call__(self, A, Y, rng=None, E=None):
tstart = time.time()
rng = np.random if rng is None else rng
A = A + rng.normal(scale=self.noise, size=A.shape) * A
X, info = self.solver(A, Y, 0, rng=rng)
info['time'] = time.time() - tstart
return self.mul_encoders(X, E), info
class _LstsqL2Solver(Solver):
"""Base class for L2-regularized least-squares solvers."""
weights = BoolParam('weights')
reg = NumberParam('reg', low=0)
solver = LeastSquaresSolverParam('solver')
def __init__(self, weights=False, reg=0.1, solver=lstsq.Cholesky()):
"""
Parameters
----------
weights : bool, optional (Default: False)
If False, solve for decoders. If True, solve for weights.
reg : float, optional (Default: 0.1)
Amount of regularization, as a fraction of the neuron activity.
solver : `.LeastSquaresSolver`, optional (Default: ``Cholesky()``)
Subsolver to use for solving the least squares problem.
Attributes
----------
reg : float
Amount of regularization, as a fraction of the neuron activity.
solver : `.LeastSquaresSolver`
Subsolver to use for solving the least squares problem.
weights : bool
If False, solve for decoders. If True, solve for weights.
"""
self.weights = weights
self.reg = reg
self.solver = solver
[docs]class LstsqL2(_LstsqL2Solver):
"""Least-squares solver with L2 regularization."""
def __call__(self, A, Y, rng=None, E=None):
tstart = time.time()
sigma = self.reg * A.max()
X, info = self.solver(A, Y, sigma, rng=rng)
info['time'] = time.time() - tstart
return self.mul_encoders(X, E), info
[docs]class LstsqL2nz(_LstsqL2Solver):
"""Least-squares solver with L2 regularization on non-zero components."""
def __call__(self, A, Y, rng=None, E=None):
tstart = time.time()
# Compute the equivalent noise standard deviation. This equals the
# base amplitude (noise_amp times the overall max activation) times
# the square-root of the fraction of non-zero components.
sigma = (self.reg * A.max()) * np.sqrt((A > 0).mean(axis=0))
# sigma == 0 means the neuron is never active, so won't be used, but
# we have to make sigma != 0 for numeric reasons.
sigma[sigma == 0] = sigma.max()
X, info = self.solver(A, Y, sigma, rng=rng)
info['time'] = time.time() - tstart
return self.mul_encoders(X, E), info
[docs]class LstsqL1(Solver):
"""Least-squares solver with L1 and L2 regularization (elastic net).
This method is well suited for creating sparse decoders or weight matrices.
"""
weights = BoolParam('weights')
l1 = NumberParam('l1', low=0)
l2 = NumberParam('l2', low=0)
def __init__(self, weights=False, l1=1e-4, l2=1e-6):
"""
.. note:: Requires `scikit-learn <http://scikit-learn.org/stable/>`_.
Parameters
----------
weights : bool, optional (Default: False)
If False, solve for decoders. If True, solve for weights.
l1 : float, optional (Default: 1e-4)
Amount of L1 regularization.
l2 : float, optional (Default: 1e-6)
Amount of L2 regularization.
Attributes
----------
l1 : float
Amount of L1 regularization.
l2 : float
Amount of L2 regularization.
weights : bool
If False, solve for decoders. If True, solve for weights.
"""
import sklearn.linear_model # noqa F401, import to check existence
assert sklearn.linear_model
self.weights = weights
self.l1 = l1
self.l2 = l2
def __call__(self, A, Y, rng=None, E=None):
import sklearn.linear_model
tstart = time.time()
Y = self.mul_encoders(Y, E, copy=True) # copy since 'fit' may modify Y
# TODO: play around with regularization constants (I just guessed).
# Do we need to scale regularization by number of neurons, to get
# same level of sparsity? esp. with weights? Currently, setting
# l1=1e-3 works well with weights when connecting 1D populations
# with 100 neurons each.
a = self.l1 * A.max() # L1 regularization
b = self.l2 * A.max()**2 # L2 regularization
alpha = a + b
l1_ratio = a / (a + b)
# --- solve least-squares A * X = Y
model = sklearn.linear_model.ElasticNet(
alpha=alpha, l1_ratio=l1_ratio, fit_intercept=False, max_iter=1000)
model.fit(A, Y)
X = model.coef_.T
X.shape = (A.shape[1], Y.shape[1]) if Y.ndim > 1 else (A.shape[1],)
t = time.time() - tstart
infos = {'rmses': rmses(A, X, Y), 'time': t}
return X, infos
[docs]class LstsqDrop(Solver):
"""Find sparser decoders/weights by dropping small values.
This solver first solves for coefficients (decoders/weights) with
L2 regularization, drops those nearest to zero, and retrains remaining.
"""
weights = BoolParam('weights')
drop = NumberParam('drop', low=0, high=1)
solver1 = SolverParam('solver1')
solver2 = SolverParam('solver2')
def __init__(self, weights=False, drop=0.25,
solver1=LstsqL2nz(reg=0.1), solver2=LstsqL2nz(reg=0.01)):
"""
Parameters
----------
weights : bool, optional (Default: False)
If False, solve for decoders. If True, solve for weights.
drop : float, optional (Default: 0.25)
Fraction of decoders or weights to set to zero.
solver1 : Solver, optional (Default: ``LstsqL2nz(reg=0.1)``)
Solver for finding the initial decoders.
solver2 : Solver, optional (Default: ``LstsqL2nz(reg=0.01)``)
Used for re-solving for the decoders after dropout.
Attributes
----------
drop : float
Fraction of decoders or weights to set to zero.
solver1 : Solver
Solver for finding the initial decoders.
solver2 : Solver
Used for re-solving for the decoders after dropout.
weights : bool
If False, solve for decoders. If True, solve for weights.
"""
self.weights = weights
self.drop = drop
self.solver1 = solver1
self.solver2 = solver2
def __call__(self, A, Y, rng=None, E=None):
tstart = time.time()
Y, m, n, d, matrix_in = format_system(A, Y)
# solve for coefficients using standard solver
X, info0 = self.solver1(A, Y, rng=rng)
X = self.mul_encoders(X, E)
# drop weights close to zero, based on `drop` ratio
Xabs = np.sort(np.abs(X.flat))
threshold = Xabs[int(np.round(self.drop * Xabs.size))]
X[np.abs(X) < threshold] = 0
# retrain nonzero weights
Y = self.mul_encoders(Y, E)
for i in range(X.shape[1]):
nonzero = X[:, i] != 0
if nonzero.sum() > 0:
X[nonzero, i], info1 = self.solver2(
A[:, nonzero], Y[:, i], rng=rng)
t = time.time() - tstart
info = {'rmses': rmses(A, X, Y), 'info0': info0, 'info1': info1,
'time': t}
return X if matrix_in else X.flatten(), info
[docs]class Nnls(Solver):
"""Non-negative least-squares solver without regularization.
Similar to `.Lstsq`, except the output values are non-negative.
"""
weights = BoolParam('weights')
def __init__(self, weights=False):
"""
.. note:: Requires
`SciPy <http://docs.scipy.org/doc/scipy/reference/>`_.
Parameters
----------
weights : bool, optional (Default: False)
If False, solve for decoders. If True, solve for weights.
Attributes
----------
weights : bool
If False, solve for decoders. If True, solve for weights.
"""
import scipy.optimize # import here too to throw error early
assert scipy.optimize
self.weights = weights
def __call__(self, A, Y, rng=None, E=None):
import scipy.optimize
tstart = time.time()
Y, m, n, d, matrix_in = format_system(A, Y)
Y = self.mul_encoders(Y, E)
X = np.zeros((n, d))
residuals = np.zeros(d)
for i in range(d):
X[:, i], residuals[i] = scipy.optimize.nnls(A, Y[:, i])
t = time.time() - tstart
info = {'rmses': rmses(A, X, Y), 'residuals': residuals, 'time': t}
return X if matrix_in else X.flatten(), info
[docs]class NnlsL2(Nnls):
"""Non-negative least-squares solver with L2 regularization.
Similar to `.LstsqL2`, except the output values are non-negative.
"""
reg = NumberParam('reg', low=0)
def __init__(self, weights=False, reg=0.1):
"""
.. note:: Requires
`SciPy <http://docs.scipy.org/doc/scipy/reference/>`_.
Parameters
----------
weights : bool, optional (Default: False)
If False, solve for decoders. If True, solve for weights.
reg : float, optional (Default: 0.1)
Amount of regularization, as a fraction of the neuron activity.
Attributes
----------
reg : float
Amount of regularization, as a fraction of the neuron activity.
weights : bool
If False, solve for decoders. If True, solve for weights.
"""
super(NnlsL2, self).__init__(weights=weights)
self.reg = reg
def _solve(self, A, Y, rng, E, sigma):
tstart = time.time()
# form Gram matrix so we can add regularization
GA = np.dot(A.T, A)
GY = np.dot(A.T, Y)
np.fill_diagonal(GA, GA.diagonal() + A.shape[0] * sigma**2)
X, info = super(NnlsL2, self).__call__(GA, GY, rng=rng, E=E)
t = time.time() - tstart
# recompute the RMSE in terms of the original matrices
info = {'rmses': rmses(A, X, Y), 'gram_info': info, 'time': t}
return X, info
def __call__(self, A, Y, rng=None, E=None):
return self._solve(A, Y, rng, E, sigma=self.reg * A.max())
[docs]class NnlsL2nz(NnlsL2):
"""Non-negative least-squares with L2 regularization on nonzero components.
Similar to `.LstsqL2nz`, except the output values are non-negative.
"""
def __call__(self, A, Y, rng=None, E=None):
sigma = (self.reg * A.max()) * np.sqrt((A > 0).mean(axis=0))
sigma[sigma == 0] = 1
return self._solve(A, Y, rng, E, sigma=sigma)