optimization 6 - neural networks corollaries

Writing up the mathematical derivation of the backpropagation got me thinking about what else can be done with it. So I've decided to write these ideas down before I forgot. Not all of them, or maybe none of them, might pan out, I haven't tried yet. What bothers me is that some of these ideas should give some substantial improvements, and they're not that hard to think of, so if they do give substantial improvements, someone would have already thought about them and they would be used all over the place. But I think they aren't. So either they're wrong and don't give the improvements, or maybe they are used and only I in my ignorance don't know about it, or maybe they really are new and worth trying.

First, I think I've been wrong thinking that doing a gradient descent step based on all the training data is too much work. If we're going to minimize the error from all the training cases, then the function to minimize will be the sum of errors from all the cases. And if we want to compute a gradient, the gradient of the sum is equal to the sum of gradients of individual parts. And we can apply the exact same principle to the whole backpropagation. So doing this compound step won't even require that much memory: just a single gradient value keeping the sum for each weight. We'd backpropagate as usual, and as we compute the gradients, we will add them to the sums. But not change the weights yet. Only after a whole pass of training is done, we'd use the computed gradients to change the weights once.

But why bother with it? Well, there are methods to select a larger gradient step rate that will get us closer to the goal in one step (such as ISTA with its Lipschitz coefficient that is used for the estimation). Just scale it down by the number of layers, to account for each layer trying to do its own equal share of adjustment. And then we can also use the momentum methods, such as FISTA (it stands for Fast ISTA) to accelerate further, for each neuron. If this works right, this should accelerate the training big time, at least for the ReLU activation function that doesn't mess up the gradients much. So if nobody uses it, there should be something wrong with this idea. It will be interesting to try. Maybe something does get lost when many small steps get replaced by one big cumulative step.

Potentially, a larger ISTA-like step rate would work on the small individual steps too, just scale the steps down more, by the number of training cases. However there is a catch in this approach: if each step takes us closer to the optimum of one training case, then the last training case will have the strongest effect, while the effect from the first training case will be most diluted. Adding them up equally, trying to do its step from the same initial position, evens them out. If the steps are small anyway, this difference doesn't matter, but if we're aiming for the larger steps, maybe it would. On the other hand, if all the cases are pushing the outputs in about the same direction anyway, doing many small separate steps would follow the gradient curve better than one large step (because, remember, we would see the zigzags with the large steps). And if we do many training cases, and repeat them in many passes, there would not be any relative advantages between the cases, they would be constantly trading spaces, the first one again becoming the last one (except on the last pass). So maybe that's why nobody accumulates the gradients. But even then I still think that a larger Lipschitz-based rate can be used, and maybe even the momentum can be carried between the passes.

The cumulative gradient for the last layer can also be used to estimate one more thing: how close we are to the optimum, and so whether it's time to stop. The squared error itself can also be used for this estimation but the trouble with it is that we don't know, how close to zero can it ever potentially get for a given configuration of a neural network. If there are many more input data than neurons (which is usually the case) then the fit won't be perfect, and we don't know in advance, how good could it be. The gradient in the last layer will be gradually dropping (more or less) as we're getting closer to the optimum, so it would be a good indicator of how far we still are. Together the gradient and the squared error can also be used to detect overfitting, when the error near optimum becomes too small. Well, maybe the cumulative gradient for the last layer, whose parts are generated not all from the same starting point but over one pass of the moving network, would work just as good.

Another problem with a possible solution is what if the random seeding of the weights generates two neurons in the same layer that are very close to each other, or even identical? They will try to stay close, and together will be a waste of space, something I've seen on my very small test examples. The solution might be that we don't have to use the same step rate for all the neurons. If we use a slightly different rate for different neurons in the layer, it would gradually push them apart, and the farther apart they get, the more should be the difference in the pressures from the gradients, so this should naturally spread them in a more useful way. But each one would still be pushed towards achieving the global optimum, just at a different speed.

This got me thinking about yet another problem, of what if a neuron with ReLU becomes dead? I.e. none of the inputs produce a positive value, and after ReLU activation, the neuron always producing 0. One possibility that I didn't notice before is that the changes in the lower-level neurons might eventually resurrect the dead neuron in a natural way by changing its inputs. Except for the neurons of the level 0, these will never see the changes of inputs. So maybe it's worth waiting longer for the higher-layer neurons to become resurrected naturally before resurrecting them forcibly by reinitialization.

And it's been bothering me that I did this reinitialization by generating new random numbers, which is both unpredictable, and as an immediate effect messes up the quality of the whole network. So maybe a way to avoid the unpredictability is by adding the weights from another neuron at a fixed relative position to this one. They've been both random to start with, so now we're adding together pairs of random numbers, which should still give random numbers. If the technique above for separating the equal neurons works well, it can even be done without disturbing the quality of the network: copy the values from another network of the same layer, and split the weight corresponding to that other neuron in each neuron of the next layer into two parts (either equal, or maybe better not) and put the other part into the weight corresponding to the resurrected neuron. Initially they together will produce the same input to the next layer as the old neighboring neuron did, and then hopefully the backpropagation pressure will be pushing them further apart.

Or yet another solution to the unpredictability might be to simply  change sings of all the weights in a dead neuron, then it will be firing on all the inputs! And to avoid disruption, set all its connected weights in the next layer to 0.

Yet another problem seemed to be that the gradient is proportional to either the input value (for the gradient of weight), or to the weight (for the input value), and so if one becomes 0, the other one stops changing. If both ever become 0, they get stuck there forever. But if at least one of them is not 0, it might be able to push the other one out, which would also provide a positive reinforcement to itself, so they would bootstrap up and it really isn't a problem. But maybe this doesn't bootstrap too well, or people would not be recommending to make the initial random values larger. So maybe a solution there would be to forcibly change a value that is too close to 0 to a larger one for the purpose of the gradient computation. It will overestimate the gradient, but that will either be a good thing, or the next pass (or even step) will restore equilibrium by pushing the values back to 0. And we can also keep some hysteresis, and pull the gradient artificially up only if it would move the point in the same direction as before.

A kind of long list of things to try, next I'll need to find time for it all!

optimization 5 - backpropagation further backwards

In the previous post we've got the last layer, how do we backpropagate from it? Backpropagation fundamentally means, how do we modify the preceding layer to change the vector v? Which can be done by doing a gradient step on v. So they turn the expression around and say that now v is the variable and x is the constant (yeah, my attempt from the last post to avoid having x as a constant was ultimately in vain, but it got us through a good deal of description). Applying the symmetrical logic we get:

df/dv[k] = 2 * x[k] * (output - b) )

Note that here x[k] has the value from before it has done its own gradient step. The estimation of the gradient by v for the next step can be improved a little by using the new value of x and recomputing the new value of output. But since the step is small, it probably doesn't matter much, and computing the new output value is extra work.

In FloatNeuralNet code this happens in the line

    wsigma += weights_[weightpos] * s_next[i];

where s_next[i] was previously set to (2 * (out - outputs[i])). The addition is there because of the next part: each neuron in the inner layers affects multiple neurons in the next layer. 

To show what is going on, let's consider a small example with a neuron of a layer-before-last with 3 inputs, that feeds into two neurons of the last layer, in beautiful ASCII art:

                                        +--> *x[1][1]----V
                                        |    *x[1][2]--> sum = f1 -> -b[1]
u[1] -> *y[1]----V                      |    *x[1][3]----^
u[2] -> *y[2]--> sum = s -> act = v[1] -+
u[3] -> *y[3]----^                      |
                                        +--> *x[2][1]----V
                                             *x[2][2]--> sum = f2 -> -b[2]
                                             *x[2][3]----^

Here I've used the variable names u and y for the input and weights of the new neuron. After summation it produces the value s, and after the activation function the value v[1]. This value v[1] is then fed as the first input of the last-layer neurons. The other inputs don't matter here, so they are not shown. The weights of the last-layer neurons are still x, but now with an extra index to tell to which neuron do they belong. The results of the last-layer neurons are f1 and f2, that are compared with the training outputs b[1] and b[2].

Let's consider for now that the activation function just passes its argument through (such as ReLU for the positive values), and so v[1] = s.

The full error (AKA loss) function f that we minimize will be

f = f1 + f2

The partial derivative by v[1] (AKA the first dimension of the gradient of f) will also be a linear combination:

df/dv[1] = df1/dv[1] + df2/dv[1] =
    = 2 * x[1][1] * (output[1] - b[1]) ) + 2 * x[2][1] * (output[2] - b[2]) )

So that's where the addition comes from, the derivatives from the last-layer neurons are added up to get the derivative of the whole function that we're minimizing.

In the inner layers there are no squares on the outputs, so the derivatives are much easier to compute. If the activation function is pass-through then:

dv[1]/dy[k] = u[k]

If the activation function is not pass-through, its derivative gets included simply as a factor:

dv[1]/dy[k] = u[k] * (dv[1]/ds)

In FloatNeuralNet code this translates to the lines:

wsigma += weights_[weightpos] * s_next[i]; // in an inner loop
s_cur[prev] = der * wsigma;
// s_cur gets moved to s_next for the next iteration on the preceding layer
Value gradient = input * s_next[i];
Value newwt = saturate(weights_[weightpos] - rate * gradient);

And then it continues propagating in the same way through more and more layers.

An interesting consequence of all this is that on each layer there are two gradient descents going on: once by adjusting the weights of this layer, once by adjusting the previous layer. So if we really wanted to compute the perfect match for the weights of this layer and did it, then right after that the adjustment to the previous layers would overshoot the minimum towards the other side. So it's actually good that the descent rate value is on the small size, and it gets amplified by the number of layers:

The last layer descends by rate, and passes the gradient to the previous layer;
The next-to-last layer descends by rate, and passes the gradient to the previous layer;
... and so on..

So the effective descent rate is really(rate*n_layers) , which means that the shallow neural nets can use a higher rate value than the deep ones. Well, I guess a counteracting force might be that on the way back the gradients might be getting more and more muddled and closer to 0, so the same rate in the earlier layers will produce a smaller effect.  Apparently, that's why in the very deep networks (over 100 layers) they make some layers double-connected, feeding both the next layer and another layer that is maybe 20 or 50 layers later, to bring in the stronger gradients. I wonder if a similar effect can be obtained by using a larger descent rate in the earlier layers.

optimization 4 - backpropagation for the last layer

As you might have guessed already from discussion of gradient descent, there is a direct relation between non-linear optimization and backpropagation. Previously I've pointed to https://www.researchgate.net/publication/266396438_A_Gentle_Introduction_to_Backpropagation as a good intro reading. Having the background optimization and looking at the neural network training problem as an optimization problem can explain, where does all this logic come from. Here is this explanation:

First, about the variable names. Traditionally, the linear regression stage in a neuron is described as:

output = sum(x[i] * w[i])

where I'll use the programming-like square-bracket indexes instead of the math-like subscript indexes; x is a vector of inputs and w is a vector of weights. By the way, what we call in programming the indexes or elements of a vector, in mathematics they call dimensions, sometimes not differentiating between the index and element, and using the word "dimension" for both. I'll use this word sometimes too. 

The thing though is that when doing the training, in this naming x is a constant and w is our variable that we're trying to optimize to achieve the training. Mathematically, it doesn't matter which letter is used for what, but I personally am used to x being a variable, so seeing it as a constant confuses me and makes me make more errors. So I'll use a different assignment of the names:

output = sum(v[i] * x[i])

where v is a vector of inputs and x is a vector of weights.

And let's ignore the non-linear activation function after the linear regression at the output side of the neuron for the moment. My experience shows that you can't use the ReLu function at the output of the last layer anyway, because the negative side won't backpropagate properly, and the Gentle Introduction also says that the activation function is not used after the last layer (and before this post is over, you'll see why). Thus let's start the consideration from the last layer that doesn't have this extra complexity.

We want to train the neuron to minimize the squared error on all the examples. I.e. if we have a set of N training examples, so that j-th example is (v[j], b[j]) where v is the vector of inputs and b is the expected output, then we want to minimize

min(f(x)),
f(x) =  sum{1 <= j <= N}( (ouput(v[j], x) - b[j])^2 

and if we expand the ouput() function and assume that the neuron has K inputs, that expression becomes:

f(x) = sum{1 <= j <= N}( (sum{1 <= i <= K}(v[j][i]*x[i]) - b[j])^2 )

If we ignore any boundaries, minimization would require just finding the values of x where the function is minimal, and since the function is quadratic and the shape of the curve will be a multidimensional parabola, there will be exactly one minimum that will be the global minimum. At this minimum the gradient (i.e. the multi-dimensional derivative) will be equal to 0, so it's even possible to find the minimum analytically: each dimension of the gradient will be linear (because the original function is quadratic, and a derivative loses one power), so equaling it to 0 will produce a system of linear equations that can be solved. 

Each dimension k of the gradient (the proper name for it is the "partial derivative of f() by x[k]") will be:

df/dx[k] = sum{1 <= j <= N}( 2 * v[j][k] * (sum{1 <= i <= K}(v[j][i]*x[i]) - b[j] )) 

where j is the index of the  training example, i is the index of the input of the neuron, and k is the index of the gradient's dimension (which has the same range 1 <= k <= K as the index of inputs).

An exact solution of a system of linear equations might be not so bad complexity-wise for a small number of inputs, but can get much worse for a large number of inputs. In that case it can be approximated by a gradient descent method, and with momentum methods it can be not so bad, say 100 steps of descent, each of quadratic time complexity, to get a decent approximation. The trouble though is that gives the perfect solution for only one neuron of the last layer, and then we'll do the backpropagation to the previous layers, and the inputs v will shift, and we'll need to do the same thing again. And also while we compute it, we'll need to remember the outputs of all the layers for all the examples to do the backpropagation. So when looking at the whole training problem, the complexity balloons very fast.

Instead people have come up with an inexact approximation. They said, let's do one training case at a time, it will give us some gradient step, the next case will give some other gradient step, and after all the cases are done, we'll get somewhere reasonably close to where we'd have been with the perfect soluton. And a small difference won't matter much because the inputs v will be shifting on the next step anyway. And we'll need to keep much less state at a time, just the outpus of all the layers for a single case, so every which way it's much cheaper. Because of these arbitrary wanderings in the middle of each training pass, they call it the stochastic gradient descent.

So OK, since we're looking at only one training case now, let's drop the index j, and the functions become:

f(x) = (sum{1 <= i <= K}(v[i]*x[i]) - b)^2 

df/dx[k] = 2 * v[k] * (sum{1 <= i <= K}(v[i]*x[i]) - b) )

I personally find the more concrete examples easier to understand and compute, so here is how I came up with the formulas above: I started with an example of a neuron with 3 inputs (i.e. K=3). Then the specific formula will be:

f(x) = (v[1]*x[i] + v[2]*x[2] + v[3]*x[3] - b)^2

To compute the gradient dimension, we'd need to expand this square, but we can note than in taking the derivative each term that doesn't contain x[k] will disappear, so we need only to bother with the terms that contain at least one x[k]. So for example for k=1 we can use the function approximation

f(x) = v[1]^2*x[1]^2 + 2*v[2]*x[2]*v[1]*x[1] + 2*v[3]*x[3]*v[1]*x[1] + (...)

df/dx[1] = 2*v[1]^2*x[1] + 2*v[2]*x[2]*v[1] + 2*v[3]*x[3]*v[1] =
= 2*v[1] * (v[1]*x[1] + v[2]*x[2] + v[3]*x[3])

And the same symmetrically for the other dimensions of x.

Alright, we've got the gradient. We can do one more simplification. Note that the sum inside it is equal to the output of the linear regression:

output = sum{1 <= i <= K}(v[i]*x[i]) 

So we can use it to compute the gradient:

df/dx[k] = 2 * v[k] * (output - b) )

After the gradient is computed, applying it is trivial. In the FloatNeuralNet code https://sourceforge.net/p/triceps/code/1707/tree/trunk/cpp/nn/FloatNeuralNet.cpp the gradient values are computed in two steps, first:

    s_cur[i] = 2 * (out - outputs[i]);

then in the next loop s_cur gets moved into s_next and

    Value gradient = input * s_next[i];
    Value newwt = saturate(weights_[weightpos] - rate * gradient);

Another thing to note is that this computation of gradients is valid only for the simple quadratic formula. If we were to insert an activation function after the linear regression, all the math will change, and I guess computing the gradients with included sigmoid is too much work.

And the third thing to note is that batching as it's normally done is not equivalent to  computing the minimum on the whole training set. Batching combines the training examples, and hence the errors produced on them, linearly. The minimization on the whole set combines the errors quadratically. Although I guess nothing really prevents us from doing a small number of training examples separately and then combining them quadratically. It would be a quadratic batching. But I'm not sure if it's worth the trouble. If nobody does it, it probably isn't.

By now this post is growing too large, so I'll talk about the further backpropagation in the next post.

optimization 3 - non-linear programming and gradient descent

Then go the "non-linear programming problems" (one more source of the acronym NLP, and with an even more meaningless use of the word "programming"). Here the non-linear function itself can have multiple minimums and maximums, plus there still is the boundary. The boundary, by the way, can now have non-linear inequalities too, and is the more difficult part of the problem. So the problem adds the part of finding the function extremums and deciding if they are in the right direction and inside the boundaries, while the part of finding the lowest (or highest) point on the boundary still remains similar to the linear programming, only is more difficult now. Now there is no easy assumption that we're always going in the same direction downhill along the boundary, so the solution is to find all the interesting points, compare them, and find the global extremum(s). And people have figured out how to find these interesting points on the boundaries, it's called the "KKT conditions" (KKT are the initials of the three authors). It's essentially looking for where the function is going down (for minimization) when nearing the boundary, and its gradient is exactly perpendicular to the boundary. The "exactly perpendicular" gives the lowest point. Although the corners where multiple boundaries intersect represent the whole range of gradients that lies between the gradients of boundary functions intersecting at this point, and the function matching any value from this range indicates a local minimum. But this all requires analytical computations of gradients and second gradients (AKA Hessians), and then analytically solving systems of inequalities on them, and none of this is easy to do.

So here enter the methods like the gradient descent. The idea there is that if the function has no local minimums, we can find the direction in which the function grows the fastest (i.e. the gradient) and go in the opposite direction to find the global minimum (if there are multiple local minimums, we would find one of them, not necessarily the best one).  

Here an interesting question is, how far should we go on each step? If we go too far, we'll overshoot the minimum. If we go not far enough, we'll have to make a lot of steps, computing the gradient on every step, which is slow.  And here is the most interesting part: even if we take the perfect step that gets us to the minimum in the direction of the gradient (which is possible to compute analytically for the quadratic functions using the Newton's method), that generally won't be the global minimum. As it turns out, the gradients in general don't give the direction to the global minimum. It gets us to some "trough" that leads towards the minimum but usually this trough is curved. And moreover, the local minimum in the direction of gradient (i.e. if the original function represents a 3-D shape, cutting in the direction of the gradient would give us a 2-D vertical slice) won't even be at the "bottom of the trough". So then if we take the second perfect step after the first perfect step, it will go at a right angle to the first one (there is a reason and a proof to this) and overshoot the "bottom of the trough" again. So the path becomes a zigzag, with smaller and smaller steps as we near the minimum (and never quite reach it).

How about we take the direction from two steps together and continue there? As it turns out, this still won't necessarily get us to the global minimum either, if the trough curves as it usually does. But we'll end up with a longer step that will be followed by a shorter perpendicular step, ad so on. We can try to correct by adding up again, which will keep the steps in the continuation direction longer again. This is the idea behind the "momentum gradient descent", where we keep the "momentum of inertia" from the previous steps, sort of imitating a ball rolling down a hilly landscape. It still will be making shorter and shorter meaningful steps as we get closer to the minimum (as the ball slows down when it nears the low spot), but at least it auto-adapts by growing the step size as necessary at the start, building up the longer steps out of the too-short fixed ones, so we can skip the complex computations of the Newton method and just err on the short side when selecting the step size. The momentum methods still work better with the basic step size that is closer to optimum (very short basic steps cause it to "circle the drain" around the optimum a lot), but still the difference is not as tremendous as with the plain gradient descent.

optimization 2 - solving the linear programming

The classic way of solving the linear programming problems is the Simplex Algorithm. It basically starts with an arbitrary corner (a tricky part is to find any starting corner) and then goes ("pivots") from there to the next neighboring corner that is lower, and so on, until it finds the lowest corner. Only it's all done as computations on vectors and matrices. It also turns out that as we get more than 2 dimensions to intersect with the plane of the function, there can be multiple corners at the same level that is not optimal, while the whole shape still is convex. So it may take the algorithm some time to sort out though these equal corners before finding one that las a connection to a corner going further down, or it may even get stuck going in circles between the same corners (but it will know that it's not at an optimum yet). For some reason they didn't talk about the obvious solution for going in circles: remember, which corners have been already tried at the current "plateau" and don't try them the second time, try the different ones (and yes, you can see, which corners are part of the plateau, they will have the key metric equal to 0).

And as it turns out, each LP problem has a related "dual" problem that it can be mechanically converted to (and it's symmetric, the second conversion, the "dual of the dual" will give the original "primal" problem). Mathematically, the dual problem is produced essentially by transposing the matrix that represents the set of inequalities, "turning it sideways", so that each variable in the primal problem is converted to an inequality in the dual one, and each inequality in the primal problem to a variable in the dual one. Logically, the relation between the problems is like this: A typical (um, "standard form" if I remember right) LP problem can also be represented as a graph of pipelines (or trade routes) of certain capacity, with the question of how to maximize the flow between two given points by using all the possible intermediate paths to the maximum of capacity - how much do we sent through each pipe? The dual problem for it would be, if we want to cut the connection between these two points completely by blocking the pipes with least possible total capacity, which pipes do we block?

Both sides of the problem are about finding the bottleneck, the primal about using it to the max, the dual about blocking it, so if you find the bottleneck in one way, it gives you the solution for the other way too. And it turns out that the dual version may be easier to solve: the major difficulty for Simplex Algorithm is finding any valid corner as a starting point, and the conversion from primal to dual provides this starting point for the dual problem. Hence the Dual Simplex Algorithm.

optimization 1 - linear programming

This year my daughter was taking a college class on optimization, so I've had a look too. No, it's not about the software optimization, it's about optimization as a math discipline. And I've learned a good deal of new things that I want to share (and also write down a summary for myself, before I forgot too much of it). The reason why you'd want to care about it is that it has an application to the neural networks, the gradient descent algorithm used for their training is one of the optimization methods, and there are interesting consequences to think about. I'll get to that part eventually (it will take multiple posts), but first things first.

They've started with the "linear programming" problems. I've been thinking "oh, I know this stuff, it's all recursion + memoization". But no, it turns out that they use this name for something completely different. In programming we call "linear programming" the discrete problems where we need to find a best-fitting combination of the discrete units. Not much programming but at least some in all that recursive trying. What they do in optimization is they call "linear programming" the analog problems where the expressions are linear.

So you have a number of variables - call them x, y, z, ..., or x_1, x_2, x_3, ..., or a vector x (this subject is ripe on vectors and matrices). And you have a function on them that you want to maximize or minimize. A linear function (that's the source of "linear" in "linear programming") would be:

  c_1 * x_1 + c_2 * c_2 + c_3 * x_3 + ...

Well, the solution would be easy, to minimize it, send all the values of x to minus infinity (or if some c < 0 then those x to plus infinity), right? But the problems also come with conditions on x, given as inequalities, which in the "linear programming" are also linear expressions. A simple example would be

  x_1 > 0

or a more complex example

  3 * x_1 - 5 * x_2 < 5

So instead of going to infinity, we'll hit the boundaries of conditions, and the problem is pretty much about which boundaries will be most restrictive.

In visual terms, the conditions form a polyhedral shape (which might be enclosed or open on some side). In a 2-dimensional form (with only x_1 and x_2) we can represent it as a 2-dimensional fence, and then the function will be some inclined plane. If we build a model of this out of plywood and let a ball roll on the plywood sheet that represents the inclined plane of the function, the ball would stop against the fence at the lowest point and solve the minimization problem (BTW, the minimization and maximization problems can be turned into each other by simply changing the signs on everything; including turning around the "less" and "greater" signs on the boundaries). There would always be a corner (or maybe a whole straight line) that is at the lowest level. Or if the fence happens to be open at the lowest point, the ball would roll out to "minus infinity", that's called an "unbounded problem". It could also be that the conditions contradict each other, and when we try to bulid the fence, we could not do it - that's a "problem with no solution". But since each condition represents an infinitely long straight line, the fenced area cut off by all these lines is guaranteed to be convex (a "convex polyhedral"), with no strangely-shaped nooks and crannies.

In the full multi-dimensional form the conditions represent the planes in that multidimensional space, that cut and split that space, and form this multidimensional "fence". This fence would still be a convex polyhedral, and the solution would be where the condition plane intersects this "fence" at the lowest point in the condition dimension.

For solving this, they have two slightly different "standard" forms of presenting these equations, one is called the "standard form" and another the "canonical form". I guess, two different people approached these problems differently, so now there are forms for both approaches, but the problems can be converted easily between these forms.

small servers

 Pretty regularly we (the programmers in general) do make the small servers that sit on some random port and execute commands from the clients. With no security. We rely on the obscurity of the randomly (but not really randomly, since we usually choose it explicitly) chosen port, and usually also on the relative privacy of the development machines for security. And ideally, firewalls that forbid connections to random ports from the outside. But what if the machine is not very private (such as the tool gets run on a production machine), and the firewall doesn't forbid the connection? Someone else might connect. Then of course there is the obscurity of the tool's protocol, but if the tool becomes popular, its protocol becomes widely known. And then if the tool can access files, it easily becomes a target for attacks. So a little obscure tool can outgrow its breeches and become a security hole.

But it seems to be not that hard to fix with a helper library that would check the password on connection. It can automatically generate the password and put it into a file, and a client would read it from the file and supply on connection (like .Xauthority does). Maybe there even is already an existing library that does this?

constants in Perl

When writing Triceps code, I've been whining that it would be great to have symbolic constants in Perl, this would allow to identify the fields of Triceps records by symbolic names that would translate at compile time into integers. Well, it turns out that there are constants in Perl and have been for more than a decade now. They've been there when I've started writing Triceps, just I wasn't aware of them. They're used like this:

use constant X => 1;
use constant Y => 2;

$a[X] = $a[Y] + 1;

They're referenced by names without any prefix, and can be scoped in packages and exported from them like everything else. Duh.

Jerry Baulier

 I don't often look at Linkedin, and even less often look at their feed (it's extremely boring), but I've opened it today and right at the top of the feed there is the news that Jerry Baulier who was the CTO of Aleri has died on April 29th. Apparently he had retired from SAS last Fall.

first attempts at training FloatNeuralNet

 After experimenting with the "intelligent design", I've tried out training the same network from a random state to do the same thing, computing the square function. Yes, it's not a typical usage, but it's a surprising example, I think automatically approximating an analog function  is quite impressive.

And it worked. Not every time, depending on the initial random state, but it when it worked, it worked quite decently. Usually worked better with more neurons. My guess is that with more neurons there is a higher chance that some of them would form a good state. Some of them end up kind of dead but if there were more to start with, it's more likely that some useful number of them will still be good.  The ReLU function I've used has this nasty habit of having a zero derivative in the lower half, and if a neuron gets there on all inputs, it becomes dead.

Which brought the question: how do we do a better initialization? Yes, I know that the bad training from a bad random initialization is normal, but I wonder, could it be done better? I've found the article: https://machinelearningmastery.com/rectified-linear-activation-function-for-deep-learning-neural-networks/.  Initialization-wise it boils down to initializing the weights between the neurons in range +-1/sqrt(num_of_inputs), and initializing the bias to 0.1. There is also an idea to consider the ReLU derivative in the negative side as a small negative value but I don't want to touch that yet. The square root probably comes from the variance of the standard distribution. The bias of 0.1 "pulls" the value up, so that even if the weights are mostly or all negative, there is still a chance of the neuron producing a positive output that would pass through ReLU. This did improve things, but only slightly, and still pretty regularly I would get a bad training.

The next thing I've noticed was that things tend to go better if I do the first few rounds (where presenting the whole set of input data is a round) of training at a higher rate. The recommendation from the article on backpropagation was to use the training rate of 1%.  I've found that if I do the first few rounds with the rate of 10%, I get fewer bad trainings. I've eventually settled down on the first 50 rounds out of 1000 using the high rate.

BTW, I've been worrying that when the weight gets initialized with a certain sign, it would always stay with the same sign, because of the ReLU having a zero derivative in the lower part. The experiment has shown that this is not the case, the weights do change signs, and even not that rarely. I still need to look at the formulas to understand, what did I miss in them.

What next? It would be good to reclaim the dead neurons and put them to a good use. Suppose we do a round of training with the whole training set. If a neuron never produces a positive value on the whole training set, it's effectively dead, and can be reinitialized by randomization. So I've done this. It does add some overhead but most of it happens infrequently, at most once per round. The only overhead to the normal training is increasing the usage counters. Observing the effects, he weird part is that it's not like "once you get all the neurons into a live state, they stay alive". Yes, the dead neurons tend to be found near the start of training but not necessarily on the consecutive rounds, and once in a while they pop up out of nowhere on a round like 300 or 500. I've added the reclamation of dead neurons after each round, and this made a big improvement. I was even able to shrink the network to the same size as I used in the "intelligent design" version and get the good results most of the time, normally better than from my handmade weights (although not as good as handmade + training). But not all the time.

Looking at the weights in the unsuccessful cases, the problem seems to be that it happens when all the neurons in this small model get initialized too similarly to each other. And then they get pushed by the training in the same direction. This seems to be the reason why the larger training rate on the early rounds made an improvement: it helps the neurons diverge. What could I do to make the neurons more different? Well, they got the fixed bias weights of 0.1. How about changing that, since I've got another way of reclaiming the dead neurons now? I've made the bias weights random, and this reduced the frequency of bad cases but didn't eliminate them.

Another thing commonly seen in the bad trainings seems to be that the weights in level 1 tend to be pushed to a low absolute value. Essentially, the NN tries to approximate the result with almost a constant, and this approximation is not a very good one. This drives the weights in the second layer up  to absolute values greater than 1, to compensate for that constness. How about if I won't let the weights get out of [-1, 1] range, will this push the lower layer to compensate there? It did, the weights there became not so small, but overall the problem of the weights driven towards 0 is not solved yet in my example. In an article I've read in IEEE Spectrum, they re-train the networks to new input data by reinitializing the neurons with low absolute weights. So maybe the same thing can be done in the normal training.

P.S. Forgot to mention, I've also changed the initialization of weights to be not too close to 0. With the upper limit on the absolute value limited by the formula with square root, I've also added the lower absolute limit at 0.1 of that.

developing FloatNeuralNet in Triceps

For a while I've wanted to experiment with the neural networks but haven't got around to it. I did some playing with them in a training course a few years ago, but there it was a black box into which you plug the training data, not the same thing as making your own.

Doing it in Perl, as I've previously done with Bayesian examples, would probably be slow. Finally I've realized that Triceps is not that much different from TensorFlow: both allow implementing the data flow machines, just with different processing elements. So nothing really stops me from adding a neural network processing element to Triceps. And it would be a good platform for experimentation, and at the same time potentially useful as production code.

My first idea was to implement the neural network as short (16-bit or even 8-bit) integers, since apparently the people say that it's good enough precision.  In reality it's not so simple. It might work for the computations in a pre-trained network but it doesn't work at all for training. In training the values can easily get out of the [-1, 1] range, and also there is a need for a decent precision: if you update the weights by 1% of the gradient, that quickly gets out of the precision that can be represented in 16 bits. On top of this there is the trouble with representation: what value do you take as 1? If INT16_MAX, it's not a power of 2, so after you multiple two numbers, you can't just shift the values to take its top bits as the result (and no, it's not the top 16 bits of an int32, it would be bits 15 to 30 of int32). Before you do a shift, you have to do another multiplication, by (INT16_MAX + 1)/INT16_MAX. If you don't do it, the values will keep shrinking with consecutive multiplications. Of if you take 0x4000 to denote 1, the shifting becomes easier but you lose one bit of precision. All in all, it's lots of pain and no gain, so after a small amount of experimentation I've abandoned this idea at least for now (the vestiges of the code are still there but likely will be removed in the future), and went with the normal floating-point arithmetic. The integer implementation is one of those ideas that make sense until you actually try to do it.

The code sits in cpp/nn. It's not connected to the Triceps row plumbing yet, just a bare neural net code that works with raw numbers.

The first part I wrote was the computation with existing weights, so I've needed some weights for testing. I've remembered a story that I've read a few years ago. It's author, a scientist, was reviewing a paper by his colleagues that contained a neural network trained to recognize the dependency between two variables. He quickly realized that had his colleagues tried the traditional methods of curve fitting, they'd see that they've trained their neural network to produce a logarithmic curve. So how about I intelligently design a curve? Something simpler than a logarithm, a parabola: y = x². It's pretty easy to do with the ReLU function to produce a segmented line: just pick a few reference points and make a neuron for each point. I've had this idea for a few years, but It's been interesting to put it into practice.

I've picked 4 points in range [0, 1] to produce 3 segments, with x values: 0, 1/3, 2/3, 1. Then computed the angle coefficients of the segments between these points (they came to 1/3, 1, and 5/3), and went to deduce the weights. The network would be 2-level: the first level would have a neuron per segment, and the second level would add them up. The second level's weights are easy: they're all 1, and bias 0. For the first neuron of level 1, the weight of the single input is equal to the angle coefficient 1/3 and the bias weight is 0. For the second neuron, we have to take the difference between the second and first angle coefficients, because they get added up on level 2, so it's 1 - 1/3 = 2/3. And the bias has to be such that this segment "kicks in" only when the value of x grows over the point's coordinate 1/3. Except that we have to work in terms of y that gets computed by weights, not x. So the bias = x * 2/3 = 1/3 * 2/3 = 2/9. Except that we also have to subtract it to make the lower values negative and thus cut off by ReLU, so bias = -2/9. And for the third point we similarly get the weight 2/3  and bias -4/9. The final formula is:

y = ReLU(ReLU(1/3 * x) + ReLU(2/3 * x - 2/9) + ReLU(2/3 * x - 4/9))

This is not the best approximation, since it all sits above the true function, but close enough:

  0.0 ->   0.000000 (true   0.000000)
  0.1 ->   0.033333 (true   0.010000)
  0.2 ->   0.066667 (true   0.040000)
  0.3 ->   0.100000 (true   0.090000)
  0.4 ->   0.177778 (true   0.160000)
  0.5 ->   0.277778 (true   0.250000)
  0.6 ->   0.377778 (true   0.360000)
  0.7 ->   0.500000 (true   0.490000)
  0.8 ->   0.666667 (true   0.640000)
  0.9 ->   0.833333 (true   0.810000)
  1.0 ->   1.000000 (true   1.000000)

I've computed the error on these values from 0 to 1 with step 0.1, and the mean error is 0.0166667, mean squared error 0.019341. Not that bad. Next I implemented the training by backpropagation, as described in https://www.researchgate.net/publication/266396438_A_Gentle_Introduction_to_Backpropagation. So what happens when we run the training over an intelligently designed weights? They improve. After running 100 rounds with the same values 0 to 1 with step 0.1, the results become:

  0.0 ->   0.000000 (true   0.000000)
  0.1 ->   0.013871 (true   0.010000)
  0.2 ->   0.047362 (true   0.040000)
  0.3 ->   0.080853 (true   0.090000)
  0.4 ->   0.156077 (true   0.160000)
  0.5 ->   0.256511 (true   0.250000)
  0.6 ->   0.356945 (true   0.360000)
  0.7 ->   0.484382 (true   0.490000)
  0.8 ->   0.651865 (true   0.640000)
  0.9 ->   0.819347 (true   0.810000)
  1.0 ->   0.986830 (true   1.000000)

The mean error drops to 0.000367414, mean squared error to 0.00770555, a nice improvement. Running the training for 1000 rounds makes the mean error very tiny but the mean squared error doesn't drop much. What happens is that the training had shifted the line segments down where they became mostly symmetrical relative to the true curve, thus the mean positive and negative errors cancelling each other out. But there is only so much precision that you can get with 3 segments, so the mean squared error can't drop much more. In case if you're interested, the formula gets shifted by training to:

ReLU(1.001182 * ReLU(0.334513 * x - 0.015636) + 1.002669 * ReLU(0.667649 * x - 0.225438)  + 1.002676 * ReLU(0.668695 * x - 0.441155) )

The changes in weights are very subtle but they reduce the errors. Which also means that at least in this particular case the 8-bit numbers would not be adequate, and even 16 bits would not be so great, they would work only in cases that have larger errors to start with.

And yes, there are the weights greater than 1! The can moved down to 1 by reducing all the weights in their input neurons but it apparently doesn't happen very well automatically in the backpropagation. It looks like the pressure for values (i.e. gradient) in the backpropagation gets distributed between the layers, and all the layers respond to it by moving in the same direction, some getting out of the [-1, 1] range. Reading on the internets, apparently it's a known problem, and the runaway growth of values can become a real issue in the large networks . The adjustment would be easy enough by rebalancing the weights between the layers but that would probably double the CPU cost of training, if not more. So for the small models, it's probably good enough to just ignore it. Well, except for one consequence: it means that the implementation of ReLU function must be unbounded at the top, passing through any values greater than 0.

Another interesting property of ReLU function is that its derivative being 0 for negative values, means that the values with the wrong sign won't be affected by training. For example, I've also created the negative side of the parabola with the same weights in negative (and for something completely different, connected them to a separate level 2 neuron), and training the positive side doesn't affect the negative side at all.

Another potential issue would apply to any activation function. The computation of the gradient for a particular input includes:

gradient = input[i] * sigma_next[j]

Which means that if the value is 0, the gradient will also be 0, and so the weights won't change. Which in particular means that a neuron that starts with all weights at 0 will produce 0 for any input, and won't ever get out of this state through backpropagation. That's one of the reasons to start with the random weights. I've found the document https://machinelearningmastery.com/rectified-linear-activation-function-for-deep-learning-neural-networks/  that talks about various tricks of initialization, and my guess is that their trick of initializing the bias weights to a small positive value is caused by this issue with zeros. But there might be other workarounds, such as messing with the activation function for gradient computation to make it never produce 0. This is the beauty of having your own code rather than a third-party library: you can experiment with all the details of it. But I haven't messed with this part yet, there are more basic things to implement first.

my another book

 I've published my book on car setup for racing:

https://www.amazon.com/dp/B09SBRGFYL

Not really a computer related subject, but anyway. Unfortunately, it's on the evil Amazon, but unfortunately there isn't a whole lot of other options at  the moment :-(

a book on Perl

 I've been browsing at a used books store, and bought the book "Higher order Perl" by Mark-Jason Dominus.  It's been published in 2003 but is still pretty amazing. It's about doing things that people normally associate with languages like Haskell, but in Perl. But it doesn't stop there. The book goes on to show an implementation of a parser infrastructure somewhat like ANTLR (but as you can imagine, with a lot less code), which is pretty mind-boggling, and then goes on to use it for a declarative drawing system that solves linear equations to determine the positions of the elements as described in a domain-oriented language.

The book can now be downloaded for free: https://hop.perl.plover.com/#free from the web site.

sentiments

 It's pretty amazing how the development of AI research is now intertwining with psychology. Not really surprising, since if the current ML models emulate the human thinking, this is what should be happening, but the amazing part is that the development of the AI research has finally reached this point.

The IEEE Computing Edge magazine tends to be rather boring but usually contains one gem per issue. In the November 2021 issue that gem was the history of Ikonas graphics processing (I didn't realize that the raster video started to appear only in mid-1970s, before then the memory costs were prohibitively expensive, and the displays did vector graphics in hardware). In the December 2021 issue that gem is an article about the human emotions "The hourglass model revisited" (https://www.sentic.net/hourglass-model-revisited.pdf, 10.1109/MIS.2020.2992799).

They define the human emotions as a combination of 6 dimensions, 4 of them that can be positive or negative (hence forming a "hourglass", thin in the middle, growing towards extremes): sensitivity, attitude, temper, introspection, and 2 neutral that can be present or not present: expectation and surprise. Additionally they split the attitude into "towards self" and towards others".

The gradations of 4 "main" emotions from positive to negative are defined there as:

Introspection: ecstasy-joy-contentment-melancholy-sadness-grief

Temper: bliss-calmness-serenity-annoyance-anger-rage

Attitude: delight-pleasantness-acceptance-dislike-disgust-loathing

Sensitivity: enthusiasm-eagerness-responsiveess-anxiety-fear-terror

And then they define the compounds, for example "hate =  anger + fear", or "pessimism = expectation + sadness". Though strangely they define despair in the same way as pessimism and not as "expectation + grief", I guess this model has quite a bit of subjectivity.

They also define the absolute strength of an emotion (in range of -1 to +1) as a sum of strengths of its components (or we could say "of its present components") divided by the count of its present (i.e. non-0) components. 

This is an improvement on always dividing by 4, but I wonder if it could be improved by just taking a sum of highest negative and highest positive components. Their formula already shifts the weights in the same direction compared to the previous one, and arguably presence of a little bit of a side emotion would not reduce the effect of the strongest one, and maybe even add to it if both emotions are in the same direction. So maybe it should be for each direction not even

  max(abs(E_i))

but

  1-Product(1-abs(E_i))

Then I've remembered where I've seen the same problem: in the Bayesian deduction! Just rescale every sentiment from [-1, 1] to [0, 1] and consider it to be a Bayesian probability of an evidence that the sentiment is positive. And then just compose them in the Bayesian deduction, in a classic formula or with formula by odds/chances as described in https://babkin-cep.blogspot.com/2017/02/a-better-explanation-of-bayes-adaboost.html. Then the sentiments that are at 0 would naturally translate to 0.5 and won't affect the computation.

retraining the neural networks

IEEE Spectrum had printed an issue on AI, and there they've described a way to teach the old neural networks the new tricks without making them forget the old tricks: by marking the more important connections as immutable when training the network on a data set for the new kind of classification.

Which looks like a generalization of the technique of "emergents" (see https://babkin-cep.blogspot.com/2017/06/emergents.html): there everything below the top 2 layers is declared immutable and the top 2 layers are re-trained for the new classification from scratch. The idea being that the first training had taught the network, what details are important, and then re-training assembles the new classifier from these details.

But the generalization goes farther: it can look at the weights of the connections, and if the weight is close to 0 (on the range of [-1, 1]), this can't be an important connection, and is a candidate to re-train, making the model learn the new details too. They say in the article that the old classifier degrades somewhat after this re-training but that's not surprising: the set of the "emergents" has changed during the re-training while the old classifier section is still expecting the old set. It would make sense to fixate the network below the top 2 layers (the "new emergents") and do another training of the old part with that.  The newly learned emergents might even be useful to improve the classification of the previously learned categories. They actually talk in the article about doing a periodic re-training with the previously trained subsets of data but not quite in the same context.

"The Practice of Parallel Programming" available as a PDF

It kind of didn't occur to me initially to announce it here, but here we go: I've made the fully edited and typeset version of my book "The Practice of Parallel Programming" available as a PDF. It's exactly the same as the printed version, with a small intro to the Online Edition 1.0 added, and formatted to a page size that is more convenient for printing at home.

Available for download here:

https://sourceforge.net/projects/tpopp/files/tpopp-text/tpopp-pdf-210901.zip

lexical analysis with regular expressions

Today I've got this simple idea of doing the lexical analysis (i.e. tokenization) in Perl with regular expressions: write all the possible tokes as alternatives in the regular expression, and use Perl's ability to embed code into the expressions to produce the token id. For example:

$s = "++abc"; # string to tokenize
$t = 0; # token ids returned here
while ($s =~ s/^\s*([+][+](?{ $t=1; })|[+](?{ $t=2; })|[a-z]+(?{ $t=3; }))//) {
  print $t, "=", $1, "\n";
}

This looks for the tokens of "+", "++", and a simple lowercase identifier. One caveat I've found is that if one alternative is a prefix of another (such as "+" and "++"), the longer alternative must go first. But otherwise it looks pretty simple, and should be efficient. I'm probably not the first one to discover it, but I've been looking for the various lexer solutions in Perl and Python and haven't seen this one yet. Of course, it would look better in the extended syntax, and with symbolic token names.

Doing the same in Python is a bit more difficult, since it doesn't allow to execute code as part of the expression. So the parsing would be just to a string, and then matching by a string. Or then looking up the token id by a dictionary (which gets a little tricky if there could be more than one non-constant lexeme, but I guess those could be sorted in the second round if needed). Something like this:

import re
tokens = { "++": 1, "+": 2, }
# regexp can be auto-generated from reverse-ordered dict keys
lex = re.compile(r'^\s*(' + r'[+][+]' + r'|[+]' + r'|[a-z]+' + r')')
# the string to parse
s = "++abc"
m = lex.match(s)
# 3 is the token id for the identifiers that are variable
t = tokens[m.group(1)] if m.group(1) in tokens else 3
# consume the token from the string
s = s[m.span()[1]:]

It doesn't do the loop but goes through all the important steps.

P.S. In Perl, the replacement could be avoided by using the global scanning instead, scan with option /g, and use \G as the anchor for the last end of string, new start of string:

while ($s =~ /\G\s*([+][+](?{ $t=1; })|[+](?{ $t=2; })|[a-z]+(?{ $t=3; }))/g)

on Rust

 I've attended a talk on the Rust programming language, and I've had a couple of realizations.

1. The new and interesting thing in Rust is  its enforcement of "borrowing" vs "consumption" in the references, and the separation of the code that strictly follows the straight-jacketed "safe" code from the "unsafe" code.  The rest of it seems to be substandard compared to C++.

I think that most real code will end up having both the "safe" and "unsafe" parts, because the purely "safe" version is conductive only to writing the COBOL-like programs. But the separation forces you to encapsulate the unsafe parts, keep them separate and small, and then the rest of the code would just rely on them instead of spreading unsafety all over the place.

But can the same concept be imported into C++? I think it can. There actually already are the syntactic means to express the borrowing and consumption in C++. Borrowing is "const &", and consumption is "&&".  So all we need is a way to ban the plain "&" and "*" in the parts of the code that we decree "safe".

2. A thing that looks annoying in Rust is its lack of classes, it has traits instead. In this respect it's similar to Go (just as annoying) and the C++ templates (in an implicit way, and with its own even worse challenges). It obviously works but  in an annoying way.

So what is annoying? That the definitions of Rust functions end up containing a list of traits for each argument: Trait1 + Trait2 + Trait3. The same applies to Go. C++ is actually less annoying because it doesn't even require you to define the traits as such, it just uses them in templates. It pays for that by very non-obvious error messages when a trait is missing, but at least you don't have to drag this list around.

But what is a class (or an interface)? An interface is a list of methods. And each method is also a trait. So when we define an interface, what we really do is write

InterfaceA = TraitA1 + TraitA2 + TraitA3

If we could do all this "trait arithmetic", all the annoyance would go away. I wonder if Rust already supports something like that, just it wasn't mentioned in the talk?

on automatic alerting

 A few years ago I've read the book "Thinking fast and slow" by Daniel Kahneman. One of the things that stuck the chord fo rme thre was the story of how the smaller counties have the cases of both higher and lower occurrence of the diseases than the larger counties. This is very much the same as what we see with the automatic alerting: when we set up an alert for, say, request processing latency, and at night there is only one request during a time period that has an anomalously high latency and triggers an alert (and we solve this with the cludges like "don't alert if there are less than 20 points").

It looked like the number of points needs to be taken in the account, with the fewer points to be treated in a more relaxed way. But I couldn't figure out how to put it together.

Recently I've been reading some books on the probablity, and it looks like the solution happens to be well-researched: we have to treat the number of points as a sample from the full larger number of points. So if we get 1000 hits per a unit of time during the day and only 50 hits during the nighttime, we can really see these 50 hits as a sample of the full traffic, with the people sleeping serving as a filter. And then instead of looking at "is the median over N" we should look at a differnt question: given the sample, what is the probability (i.e. confidence range) that the median of the full set of data would have been over N? That should solve the issue with the false alerts pretty nicely in an automatic way. I wonder, why nobody I've heard of had laready done this?

time for more analog computers?

In case if you haven't heard about them, the analog computers were unlike the digital computers. They were the electrical models assmebled to simulate the other real-world events. For example, if you wanted to model a car suspension, you'd have an oscillator circuit to simulate a spring, and a resistor to simulate a shock absorber. They were purpose-assembled for every specific application.

Well, now we have another thing that purports to simulate the real-world systems: the units in the neural networks that are supposed to be like the real neurons. The race to cram more "neurons" into the neural networks has led to the values in the "neurons" to be represented as very short 8-bit fixed-point numbers. And the models still work. Which means that the "neurons" are quite well at tolerating the small imprecisions in the data, which is the typical issue with the analog signals. But in return the analog signals are much more resistant to large disruptions: i.e. if you change one single highest bit in a digital representation of a number, this small disruption changes the value of a number by half its whole range, something that doesn't happen with the analog signals. There already are ideas of improving the reliability of computations by making them more resistant to large errors by accepting the small errors (I think I've read about them in IEEE's Spectrum).

So, the next logical step: why not make each "neuron" in a neural network an analog machine? Instead of 8 bits, have one analog signal, and do all the multiplication and addition on the analog signals. It looks like the perfect match, a combination of these two ideas. Has anyone already done that? If not, why not? Is there a way to do that? 

P.S. It turns out, people are already doing this. I've recently read in Spectrum that people are using photonic computations for the machine learning, and that the photonic computations are limited to this niche because of their limited precision.  I can't really think of a reason for the limited precision other than these computations being analog.

P.P.S. It's 2022 now, slightly more than a year after the original post, and IEEE Spectrum is reporting that IBM is working on a way to use the memory (like Flash, but other similar types too) as an analog summator, doing the "neuron" processing in a memory cell.