Backpropagation on a matrix of functions As nicely described in this article, backpropagation for multi-layer perceptrons defines the error term for a neuron in terms of the partial derivative of the weights.

It's traditional to represent the weights between neutrons in matrix form: $w_{ij}$ representing the weight between neuroma $i$ and $j$.
Since a weight can be considered as being as constant multiplier of any input, so a possible generalisation is instead to replace the $w_{ij}$ with some univariate function $f_{ij}$.

Since. a 3 layer MLP is already a universal approximator for continuous functions, this generalisation won't be more expressive, but the hypothesis is that it might converge faster for some problems.

My question then concerns the nature of the learning algorithm for this extended formulation:

How might it be possible to define backpropagation if the weight matrix is replaced by a matrix of functions? The functions can be assumed to be differentiable as required.
The problem isn't well-defined as you've stated it.  Backpropagation involves updating the weights, i.e., finding weights that cause the neural network to behave well.  If you have no weights, then there's nothing to update.

For these reasons, the standard approach to building neural networks involves both weights and univariate functions.  In the standard perceptron architecture, there are weights $w_{ij}$ along each edge (from neuron $i$ to neuron $j$), and a univariate function $f_j$ applied to the output of each neuron (neuron $j$).  So standard artificial neural networks already incorporate the univariate function.  The univariate function $f_j$ is called an "activation function".

Standard backpropagation does work with arbitrary activation functions.  You don't need to modify or generalize it: it already deals with them.  In particular, backprop can be used with an arbitrary activation function, as long as you know its derivative.  So, I suggest you take another look at how standard backprop works.

Backprop is basically just (a) computation of the partial derivatives (the partial derivative of some objective function, with respect to the weights), computed using the chain rule, plus (b) gradient descent.  Introducing arbitrary univariate functions affects only part (a), and the chain rule from calculus tells you how to calculate the partial derivatives.