# Reservoir Computing¶

## Introduction¶

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.

## Example¶

"""
"""
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,
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,
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,
w_bias          = None,
w               = HDPy.ring_of_neurons,
w_in            = HDPy.sparse_w_in(1.5, 50, rnd_fu=np.random.normal),
)

input_dim       = reservoir_size,
output_dim      = 1,
with_bias       = True,
lambda_         = 1.0,
)

# 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)

# Initialize the reservoirs
# Propagate data through the reservoirs, no training
for res in reservoirs:
res(src[:washout])

# 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


## Reference¶

class HDPy.ReservoirNode(input_dim=None, output_dim=None, spectral_radius=0.9, nonlin_func=<ufunc 'tanh'>, reset_states=False, w_bias=None, w=None, w_in=None, **kwargs)

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).

input_dim
Number of inputs.
output_dim
Reservoir size.
Spectral radius of the reservoir. The spectral radius of a matrix is its largest absolute eigenvalue.
nonline_func
Reservoir node function, often nonlinear (tanh or sigmoid). Default is tanh.
reset_states
Reset the states to the initial one after a call to execute(). The default is False, which is appropriate in most situations where learning is online.
w

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).

w_in

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.

w_bias
Reservoir bias initialization routine. The callback is executed on one argument, the reservoir size. The default is dense_w_bias(), with the scaling 0.0 (i.e. no bias) or bias_scaling (if provided).

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.

fan_in_i
Density of the input matrix, in percent.
fan_in_w
Density of the reservoir, in percent.
input_scaling
Scaling of the input weights.
bias_scaling
Scaling of the bias.
execute(x, simulate=False)

Executes simulation with input vector x.

x
Input samples. Variables in columns, observations in rows, i.e. if N,M = x.shape then M == input_dim.
simulate
If True, the state won’t be updated.
copy()

Return a deep copy of the node.

input_dim

Input dimensions

output_dim

Output dimensions

reset()

Reset the reservoir states to the initial value.

save(filename, protocol=-1)

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.

_post_update_hook(states, input_, timestep)

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.

__call__(x, *args, **kwargs)

Calling an instance of Node is equivalent to calling its execute method.

class HDPy.PlainRLS(input_dim, output_dim, with_bias=True, lambda_=1.0)

Compute online least-square, multivariate linear regression.

The interface is based on [MDP] and [Oger] with two major differences:

1. There’s no inheritance from any [MDP] or [Oger] class
2. The error as well as a target may be passed to the train() function

The implementation follows the algorithm given in [FB98] (least-squares filter).

with_bias
Add a constant to the linear equation. This internally increases the input dimension by one, whereas the added input is always one. Generally, the target function is $$f(x) = \sum_{i=0}^N a_i x_i$$. With with_bias=True it will be $$f(x) = \sum_{i=0}^N a_i x_i + a_{N+1}$$ instead. (default True)
input_dim
Dimension of the input (number of observations per sample)
output_dim
Dimension of the output (Regression order)
lambda_

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(sample, trg=None, err=None, d=None, e=None)

Train the regression on one or more samples.

sample
Input samples. Array of size (K, input_dim)
train
Sample target. Array of size (K, output_dim)
err
Sample error terms. Array of size (K, output_dim)
__call__(x)

Evaluate the linear approximation on some point x.

save(pth)

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')
>>> f.close()
>>> (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

stop_training()

Disable filter adaption for future calls to train.

copy()

Return a deep copy of the node.

class HDPy.StabilizedRLS(input_dim, output_dim, with_bias=True, lambda_=1.0)

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(sample, trg=None, err=None, d=None, e=None)

Train the regression on one or more samples.

sample
Input samples. Array of size (K, input_dim)
trg
Sample target. Array of size (K, output_dim)
err
Sample error terms. Array of size (K, output_dim)
HDPy.sparse_reservoir(density, rnd_fu=None)

Reservoir with sparse connected neurons and random connection weights. The non-zero weights are uniformly distributed, the values drawn from a random distribution.

density
Reservoir density, given in percent (0-100).
rnd_fu
Random number generator for the weight values. Must take the scalar argument size that determines the length of the generated random number list. The default is the normal distribution.
HDPy.dense_w_in(scaling, rnd_fu=None)

Dense input matrix.

scaling
Input scaling.
rnd_fu
Callback for random number generation. Expects the argument size. Default is a uniform distribution in [-1,1].
HDPy.sparse_w_in(scaling, density, rnd_fu=None)

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.

scaling
Scaling of input weights.
density
Percentage of non-trivial connection weights, [0-100].
rnd_fu
Random number generator for the weight values. Must take the scalar argument size that determines the length of the generated random number list. The default is the uniform distribution between -1.0 and 1.0.
HDPy.dense_w_bias(scaling, rnd_fu=None)

Dense reservoir node bias. The weights are picked from a random number generator rnd_fu and scaled according to scaling.

scaling
Scale of the bias weights
rnd_fu
Random number generator for the weights. Default is normal distribution in [-1.0, 1.0].
HDPy.orthogonal_reservoir(density)

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.

density
Reservoir density, given in percent (0-100).

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.

out_size
Reservoir matrix dimension.

Reservoir setup where the nodes are arranged in a ring. Each node has exactly one predecessor and one successor.

out_size
Reservoir matrix dimension.
HDPy.reservoir_memory(reservoir, max_settling_time=10000)

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.

HDPy.find_radius_for_mc(reservoir, num_steps, tol=1.0, max_settling_time=10000, tol_settling=0.5, num_iter=100)

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.