Contents
As Reservoir Computing is uses in this module, an implementation of an Echo State Network is provided through a couple of functions. The reservoir is given through ReservoirNode. Note that although the dependency has been removed, most part of the reservoir code is copied from [Oger] and [MDP]. The reservoir has then be combined with a readout, usually a linear regression. For the offline case, a ridge regression would probably be incorporated (available in Oger.nodes.RidgeRegressionNode). Since for this module, only the online case is relevant, a recursive least-squares filter is implemented in StabilizedRLS. The reservoir and readout are easily combined by feeding the reservoir state as sample into the RLS (together with a training target).
In contrast to [Oger], the reservoir is set up through three external functions. These define how the weight matrices are initialized. The default approach is to use dense_w_in(), sparse_w_in() and dense_w_bias(). However, other reservoir setup functions are provided. Note the documentation in ReservoirNode if a different initialization has to be implemented.
To test reservoir memory, two functions are provided, reservoir_memory() and find_radius_for_mc(). Both of them measure the settling time of the impulse response, which is related to the spectral radius.
"""
"""
import HDPy
import numpy as np
import pylab
## PARAMS ##
washout = 200
num_train = 5000
num_test = 1000
reservoir_size = 100
## INITIALIZATION ##
# Reservoir
reservoir_sparse = HDPy.ReservoirNode(
input_dim = 1,
output_dim = reservoir_size,
spectral_radius = 0.9,
w_bias = None,
w = HDPy.sparse_reservoir(20),
w_in = HDPy.sparse_w_in(1.5, 50, rnd_fu=np.random.normal),
)
reservoir_orthogonal = HDPy.ReservoirNode(
input_dim = 1,
output_dim = reservoir_size,
spectral_radius = 0.9,
w_bias = None,
w = HDPy.orthogonal_reservoir(20.0),
w_in = HDPy.sparse_w_in(1.5, 50, rnd_fu=np.random.normal),
)
reservoir_ring = HDPy.ReservoirNode(
input_dim = 1,
output_dim = reservoir_size,
spectral_radius = 0.9,
w_bias = None,
w = HDPy.ring_of_neurons,
w_in = HDPy.sparse_w_in(1.5, 50, rnd_fu=np.random.normal),
)
# Readout
readout = HDPy.StabilizedRLS(
input_dim = reservoir_size,
output_dim = 1,
with_bias = True,
lambda_ = 1.0,
)
readout_orthogonal = readout.copy()
readout_ring = readout.copy()
readout_sparse = readout.copy()
# Data
def narma30(num_samples=1000):
"""30th order NARMA dataset. Copied from [Oger]_."""
system_order = 30
inputs = np.random.rand(num_samples, 1) * 0.5
outputs = np.zeros((num_samples, 1))
for k in range(system_order-1, num_samples-1):
outputs[k + 1] = 0.2 * outputs[k] + 0.04 * \
outputs[k] * np.sum(outputs[k - (system_order-1):k+1]) + \
1.5 * inputs[k - 29] * inputs[k] + 0.001
return inputs, outputs
src, trg = narma30(washout + num_train + num_test)
## TRAINING ##
setups = ('Sparse', 'Orthogonal', 'Ring of Neurons')
reservoirs = (reservoir_sparse, reservoir_orthogonal, reservoir_ring)
readouts = (readout_sparse, readout_orthogonal, readout_ring)
# Initialize the reservoirs
# Propagate data through the reservoirs, no training
for res in reservoirs:
res(src[:washout])
# Train the readout
# Propagate data through reservoir, train the readout online
for res, out in zip(reservoirs, readouts):
r_state = res(src[washout:num_train])
out.train(r_state, trg[washout:num_train])
# Test the networks
signals = []
for res, out in zip(reservoirs, readouts):
r_state = res(src[washout+num_train:])
pred = out(r_state)
signals.append(pred)
## PLOTTING ##
# Error measurement
mse = lambda sig_pred, sig_trg: ((sig_pred - sig_trg)**2).mean()
rmse = lambda sig_pred, sig_trg: np.sqrt(mse(sig_pred, sig_trg))
nrmse = lambda sig_pred, sig_trg: rmse(sig_pred, sig_trg) / sig_trg.std()
# Output and reservoir output plotting
pretty_str = "{0:<" + str(max(map(len, setups))) + "}\t{1:0.6f}\t{2:0.6f}"
print "Reservoir type\tMSE\t\tNRMSE"
for sig, lbl in zip(signals, setups):
pylab.plot(sig, label=lbl)
err_mse = mse(sig, trg[washout + num_train:])
err_nrmse = nrmse(sig, trg[washout + num_train:])
print pretty_str.format(lbl, err_mse, err_nrmse)
# Target plotting
pylab.plot(trg[washout+num_train:], 'c', label='Target')
# Show the plot
pylab.axis((0.0, 70.0, 0.0, 0.45))
pylab.legend(loc=0)
pylab.show(block=False)
>>> Sparse 0.004785 0.744121
>>> Orthogonal 0.004858 0.749770
>>> Ring of Neurons 0.004827 0.747397
Reservoir of neurons.
This class is mainly copied from [Oger], with some modifications. Mainly, the initialization is done differently (yet, compatibility maintained through additional keyword arguments) and the reservoir state can be computed without affecting the class instance (execute()‘s simulate argument).
Reservoir initialization routine. The callback is executed on two arguments, the number of nodes and the spectral radius.
The default is sparse_reservoir(), with normally distributed weights and density of 20% or fan_in_w (see below).
Reservoir input weights initialization routine. The callback is executed on two arguments, the reservoir and input size.
The default is a dense matrix dense_w_in(), with input scaling of 1.0. If fan_in_i is provided, a sparse matrix is created, instead.
For compatibility with [Oger], some additional keyword arguments may be provided. These affect the matrix initialization routines above and are only used if the respective initialization callback is not specified.
Executes simulation with input vector x.
Return a deep copy of the node.
Input dimensions
Output dimensions
Reset the reservoir states to the initial value.
Save a pickled serialization of the node to filename. If filename is None, return a string.
Note: the pickled Node is not guaranteed to be forwards or backwards compatible.
Hook which gets executed after the state update equation for every timestep. Do not use this to change the state of the reservoir (e.g. to train internal weights) if you want to use parallellization - use the TrainableReservoirNode in that case.
Calling an instance of Node is equivalent to calling its execute method.
Compute online least-square, multivariate linear regression.
The interface is based on [MDP] and [Oger] with two major differences:
The implementation follows the algorithm given in [FB98] (least-squares filter).
Forgetting factor. Controlls how much the training focuses on recent samples (see [FB98] for details). (default 1.0, corresponds to offline learning)
“Roughly speaking \(\frac{1}{1-\lambda}\) is a measure of the memory of the algorithm. The case of \(\lambda = 1\) corresponds to infinite memory.” [FB98]
Train the regression on one or more samples.
Evaluate the linear approximation on some point x.
Save the regression state in a file.
To load the state, use pickle, as show below.
>>> import numpy as np
>>> import tempfile, os, pickle
>>> fh,pth = tempfile.mkstemp()
>>> r1 = PlainRLS(input_dim=5, output_dim=2, lambda_=0.9, with_bias=True)
>>> r1.train(np.random.normal(size=(100,5)), e=np.random.normal(size=(100,2)))
>>> r1.save(pth)
>>> f = open(pth,'r')
>>> r2 = pickle.load(f)
>>> f.close()
>>> os.unlink(pth)
>>> (r2.beta == r1.beta).all()
True
>>> (r2._psi_inv == r1._psi_inv).all()
True
>>> r2.input_dim == r1.input_dim
True
>>> r2.output_dim == r1.output_dim
True
>>> r2.lambda_ == r1.lambda_
True
>>> r2.with_bias == r1.with_bias
True
Disable filter adaption for future calls to train.
Return a deep copy of the node.
Bases: HDPy.rc.PlainRLS
Compute online least-square, multivariate linear regression.
Identical to PlainRLS, except that the internal matrices are computed on the lower triangular part. This ensures symmetrical matrices even with floating point operations. If unsure, use this implementation instead of PlainRLS.
Train the regression on one or more samples.
Reservoir with sparse connected neurons and random connection weights. The non-zero weights are uniformly distributed, the values drawn from a random distribution.
Dense input matrix.
Sparse, random reservoir input. The density is exactly as defined by density, given in percent (0-100). The values are drawn from an uniform distribution.
Dense reservoir node bias. The weights are picked from a random number generator rnd_fu and scaled according to scaling.
Orthogonal reservoir construction algorithm.
The algorithm starts off with a diagonal matrix, then multiplies it with random orthogonal matrices (meaning that the product will again be orthogonal). For details, see [TS12].
All absolute eigenvalues of the resulting matrix will be identical to specrad.
Reservoir setup where the nodes are arranged in a chain. Except for the first and last node, each one has exactly one predecessor and one successor.
Reservoir setup where the nodes are arranged in a ring. Each node has exactly one predecessor and one successor.
Measure the memory capacity of a reservoir. Make sure, the reservoir is initialized. The settling time is measured, which is the number of steps it takes until the reservoir has converged after some non-zero input.
Find a spectral radius for reservoir such that the impulse response time is num_steps with tolerance tol.
This implementation linearizes the memory capacity function locally and uses binary search to approach the target value. If you look for a single value (not the whole characteristics), this function is usually faster than query_reservoir_memory().
max_settling_time sets the maximum time after which the reservoir should have settled.
tol_settling and num_iter are search abortion criteria which break, if the settling time doesn’t change too much or too many iterations have been run. If any of the criteria is met, an Exception is thrown.