Warrior_EA/AI/Network.cl
AnimateDread 4f28165cd3 fix: remove broken DFA optimizer, use plain gradient descent
The DFA (Direct Feedback Alignment) option was never a correct implementation:
it deterministically flipped the sign of half of all gradients based on
connection index parity, causing permanent gradient ascent for those weights
and guaranteed divergence. The backward pass was also incompatible with the
OpenCL/DirectML neuron model (layer.Total() == 1). This change removes all DFA
logic, including the enum value and `DfaFeedbackSignal` method, and replaces it
with plain gradient descent in all momentum update kernels. The `optimizer`
kernel argument is retained for binary compatibility but is no longer used.
2026-07-29 00:03:54 -04:00

657 lines
34 KiB
Common Lisp

// fp32 throughout - this project targets consumer GPUs (e.g. AMD Polaris/RX 580) where double
// precision runs at a small fraction of float throughput. No kernel below needs fp64 range or
// precision: weight/delta clamps stay within +-100, gradients/activations are float32-safe, and
// the Adam epsilon guard (vt > 0 ? vt : l*10) isn't a tiny-epsilon comparison sensitive to
// fp32 vs fp64 underflow. Removing cl_khr_fp64 means any accidental double slipping back in fails
// kernel compilation loudly instead of silently running at 1/16th speed.
// Tightened from 1.0e6 - see AI\Network.mqh's matching MAX_WEIGHT comment for the full rationale
// (that ceiling was so loose it never engaged before training had already gone unstable). 100.0
// matches the equivalent clamp in Dmitriy Gizlyk's reference NeuroNet.cl engine.
#define MAX_WEIGHT 100.0f
// Floor on the magnitude of a saturated tanh/sigmoid unit's activation derivative - see AI\Network.mqh's
// matching MIN_ACTIVATION_DERIVATIVE comment for the full rationale (a unit pinned near its activation
// extremes must still receive some corrective gradient, however small, rather than exactly zero).
// 2026-07-27: Increased from 1.0e-4f to 1.0e-3f to strengthen the escape signal through saturated hidden
// neurons. On the fp32 OpenCL backend, a saturated sigmoid output layer (the "neutral collapse from era 0"
// bug) propagated only this floor value backward through every hidden layer, which at 1e-4 was too weak
// to ever pull weights out of saturation. The output-layer saturation root cause is now fixed (NONE/logit
// activation instead of SIGMOID see ExpertSignalAIBase.mqh's BuildFreshTopology), but this increase
// provides a 10× stronger safety net against any remaining fp32-specific saturation in hidden layers.
#define MIN_ACTIVATION_DERIVATIVE 1.0e-3f
// Decoupled (AdamW-style) weight decay applied inside every Adam kernel below - see WarriorCPU.cpp's
// matching WEIGHT_DECAY comment for the full rationale (unbounded slow weight growth across hundreds
// of oversampled training eras causing periodic collapse-then-recover cycles). 0.01 is the standard
// AdamW default.
#define WEIGHT_DECAY 0.001f
// Per-step update clip - see WarriorCPU.cpp's matching MAX_WEIGHT_DELTA comment for the full
// rationale (weight decay alone didn't stop the collapse cycles - they turned out to be sudden Adam
// overshoot events, most likely from 5x back-to-back oversampling replay building artificially
// correlated momentum). Applied to the raw delta BEFORE it's added to the weight, unlike MAX_WEIGHT.
#define MAX_WEIGHT_DELTA 0.1f
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
__kernel void FeedForward(__global float *matrix_w,
__global float *matrix_i,
__global float *matrix_o,
int inputs, int activation)
{
int i = get_global_id(0);
float sum = 0.0f;
float4 inp, weight;
int shift = (inputs + 1) * i;
for(int k = 0; k <= inputs; k = k + 4)
{
switch(inputs - k)
{
case 0:
inp = (float4)(1, 0, 0, 0);
weight = (float4)(matrix_w[shift + k], 0, 0, 0);
break;
case 1:
inp = (float4)(matrix_i[k], 1, 0, 0);
weight = (float4)(matrix_w[shift + k], matrix_w[shift + k + 1], 0, 0);
break;
case 2:
inp = (float4)(matrix_i[k], matrix_i[k + 1], 1, 0);
weight = (float4)(matrix_w[shift + k], matrix_w[shift + k + 1], matrix_w[shift + k + 2], 0);
break;
case 3:
inp = (float4)(matrix_i[k], matrix_i[k + 1], matrix_i[k + 2], 1);
weight = (float4)(matrix_w[shift + k], matrix_w[shift + k + 1], matrix_w[shift + k + 2], matrix_w[shift + k + 3]);
break;
default:
inp = (float4)(matrix_i[k], matrix_i[k + 1], matrix_i[k + 2], matrix_i[k + 3]);
weight = (float4)(matrix_w[shift + k], matrix_w[shift + k + 1], matrix_w[shift + k + 2], matrix_w[shift + k + 3]);
break;
}
sum += dot(inp, weight);
}
switch(activation)
{
case 0:
sum = tanh(sum);
break;
case 1:
sum = 1 / (1 + exp(-clamp(sum, -50.0f, 50.0f)));
break;
case 2: // PReLU, param=0.01 - matches FeedForwardConv's case 2 below
if(sum < 0)
sum *= 0.01f;
break;
}
matrix_o[i] = sum;
barrier(CLK_GLOBAL_MEM_FENCE);
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
__kernel void CaclOutputGradient(__global float *matrix_t,
__global float *matrix_o,
__global float *matrix_ig,
int activation)
{
int i = get_global_id(0);
float temp = 0;
float out = matrix_o[i];
switch(activation)
{
case 0:
// Deliberately NOT multiplied by the tanh derivative (1-out^2): that
// factor vanishes as out approaches +-1, which is exactly where a
// TANH output neuron needs to converge for a +-1 target (e.g. the
// buy/sell extremes of a single-neuron regression head), stalling
// training right when it matters most. Using the raw (target-out)
// error here is the same fix as pairing sigmoid with cross-entropy -
// a bounded activation whose output error isn't re-damped by its
// own saturating derivative.
temp = clamp(matrix_t[i], -1.0f, 1.0f) - out;
break;
case 1:
// Also deliberately NOT multiplied by the sigmoid derivative out*(1-out) - same
// reasoning as case 0 above. This is the classification output layer (3 neurons,
// one-hot 0/1 targets - see BuildFreshTopology()); (target-out) unscaled is exactly
// the sigmoid+binary-cross-entropy gradient (the loss derivative and the sigmoid
// derivative algebraically cancel), whereas re-multiplying by out*(1-out) here is
// the MSE-with-sigmoid formula, which damps toward zero as out approaches 0 or 1 -
// exactly where a 0/1-target neuron needs to converge. With this damping (plus the
// 0.00000001 underflow guard at the saturated extremes), the classification output
// layer's effective learning rate was throttled to a small fraction of the
// regression path's, observed in practice as all 3 output neurons converging to
// (and never escaping) an identical value - softmax of three equal logits is always
// exactly 1/3 each, with ties resolved to index 0 ("Buy" wins by default).
temp = clamp(matrix_t[i], 0.0f, 1.0f) - out;
break;
case 2:
// PReLU, param=0.01 - unlike tanh/sigmoid, PReLU is unbounded and never saturates,
// so none of the clamp-to-range/derivative-skip workarounds above are needed here:
// plain (target-out) scaled by the derivative at this output is exact.
temp = (matrix_t[i] - out) * (out >= 0 ? 1.0f : 0.01f);
break;
default:
// NONE (raw logits - see the 3-class softmax classification output layer in
// BuildFreshTopology()/NativeActivationCode()'s NONE case, which maps here): identity
// activation, derivative 1 everywhere, so the plain (target-out) error passes through
// unscaled - same formula CNeuron::calcOutputGradients() already uses unconditionally
// on the un-accelerated CPU fallback path. Leaving this case unhandled left `temp` at
// its 0.0 initializer, i.e. every softmax-classification output neuron got a
// permanently zero gradient on this backend - which, since calcHiddenGradients()
// propagates backward FROM the output layer's gradient, silently froze the entire
// network (not just the output layer) at its random initial weights whenever an
// OpenCL-accelerated tier was active.
temp = matrix_t[i] - out;
break;
}
matrix_ig[i] = temp;
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
__kernel void CaclHiddenGradient(__global float *matrix_w,
__global float *matrix_g,
__global float *matrix_o,
__global float *matrix_ig,
int outputs, int activation)
{
int i = get_global_id(0);
float sum = 0;
float out = matrix_o[i];
float4 grad, weight;
int shift = (outputs + 1) * i;
for(int k = 0; k < outputs; k += 4)
{
switch(outputs - k)
{
case 0:
grad = (float4)(matrix_g[k], 0, 0, 0);
weight = (float4)(matrix_w[shift + k], 0, 0, 0);
break;
case 1:
grad = (float4)(matrix_g[k], matrix_g[k + 1], 0, 0);
weight = (float4)(matrix_w[shift + k], matrix_w[shift + k + 1], 0, 0);
break;
case 2:
grad = (float4)(matrix_g[k], matrix_g[k + 1], matrix_g[k + 2], 0);
weight = (float4)(matrix_w[shift + k], matrix_w[shift + k + 1], matrix_w[shift + k + 2], 0);
break;
default:
grad = (float4)(matrix_g[k], matrix_g[k + 1], matrix_g[k + 2], matrix_g[k + 3]);
weight = (float4)(matrix_w[shift + k], matrix_w[shift + k + 1], matrix_w[shift + k + 2], matrix_w[shift + k + 3]);
break;
}
sum += dot(grad, weight);
}
switch(activation)
{
case 0:
sum = clamp(sum + out, -1.0f, 1.0f) - out;
sum = sum * fmax(MIN_ACTIVATION_DERIVATIVE, 1 - pow(out, 2.0f));
break;
case 1:
sum = clamp(sum + out, 0.0f, 1.0f) - out;
sum = sum * fmax(MIN_ACTIVATION_DERIVATIVE, out * (1 - out));
break;
case 2:
// PReLU, param=0.01 - unbounded, non-saturating, so the plain chain rule applies
// directly with no clamp/implied-target reformulation needed (contrast with the
// tanh/sigmoid cases above, which exist specifically to dodge saturation).
sum = sum * (out >= 0 ? 1.0f : 0.01f);
break;
}
matrix_ig[i] = sum;
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
__kernel void UpdateWeightsMomentum(__global float *matrix_w,
__global float *matrix_g,
__global float *matrix_i,
__global float *matrix_dw,
int inputs, float learning_rates, float momentum, int optimizer)
{
int i = get_global_id(0);
int j = get_global_id(1);
int wi = i * (inputs + 1) + j;
// `optimizer` is retained (unused) only to keep this kernel's argument list binary-compatible with
// the already-deployed WarriorDML.dll/WarriorCPU.dll exports. It used to select an index-parity sign
// flip on the gradient ("DFA"), which put half of every weight tensor into permanent gradient ASCENT
// - see ENUM_OPTIMIZATION's comment in AI\Network.mqh. Always plain descent now; callers pass 0.
float delta = learning_rates * (matrix_g[i] * (j < inputs ? matrix_i[j] : 1)) + momentum * matrix_dw[wi];
matrix_dw[wi] = delta;
// Unlike UpdateWeightsConvMomentum/UpdateWeightsConvAdam/LSTM_UpdateWeightsAdam below, this dense
// kernel (and UpdateWeightsAdam right after it) had never gotten the MAX_WEIGHT clamp - a gradient
// spike could drive a weight to +-Infinity with nothing to stop it, and the next Adam step turns
// Infinity into NaN (Inf/Inf), which then propagates through every FeedForward sum touching that
// weight and never recovers. Dense layers are exactly where the classification output head lives,
// so this is the most consequential of the three missing clamps, not the least.
matrix_w[wi] = clamp(matrix_w[wi] + delta, -MAX_WEIGHT, MAX_WEIGHT);
};
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
__kernel void UpdateWeightsAdam(__global float *matrix_w,
__global const float *matrix_g,
__global const float *matrix_i,
__global float *matrix_m,
__global float *matrix_v,
const int inputs, const float l, const float b1, const float b2)
{
const int i = get_global_id(0);
const int j = get_global_id(1);
const int wi = i * (inputs + 1) + j * 4;
float4 m, v, weight, inp;
switch(inputs - j * 4)
{
case 0:
inp = (float4)(1, 0, 0, 0);
weight = (float4)(matrix_w[wi], 0, 0, 0);
m = (float4)(matrix_m[wi], 0, 0, 0);
v = (float4)(matrix_v[wi], 0, 0, 0);
break;
case 1:
inp = (float4)(matrix_i[j], 1, 0, 0);
weight = (float4)(matrix_w[wi], matrix_w[wi + 1], 0, 0);
m = (float4)(matrix_m[wi], matrix_m[wi + 1], 0, 0);
v = (float4)(matrix_v[wi], matrix_v[wi + 1], 0, 0);
break;
case 2:
inp = (float4)(matrix_i[j], matrix_i[j + 1], 1, 0);
weight = (float4)(matrix_w[wi], matrix_w[wi + 1], matrix_w[wi + 2], 0);
m = (float4)(matrix_m[wi], matrix_m[wi + 1], matrix_m[wi + 2], 0);
v = (float4)(matrix_v[wi], matrix_v[wi + 1], matrix_v[wi + 2], 0);
break;
case 3:
inp = (float4)(matrix_i[j], matrix_i[j + 1], matrix_i[j + 2], 1);
weight = (float4)(matrix_w[wi], matrix_w[wi + 1], matrix_w[wi + 2], matrix_w[wi + 3]);
m = (float4)(matrix_m[wi], matrix_m[wi + 1], matrix_m[wi + 2], matrix_m[wi + 3]);
v = (float4)(matrix_v[wi], matrix_v[wi + 1], matrix_v[wi + 2], matrix_v[wi + 3]);
break;
default:
inp = (float4)(matrix_i[j], matrix_i[j + 1], matrix_i[j + 2], matrix_i[j + 3]);
weight = (float4)(matrix_w[wi], matrix_w[wi + 1], matrix_w[wi + 2], matrix_w[wi + 3]);
m = (float4)(matrix_m[wi], matrix_m[wi + 1], matrix_m[wi + 2], matrix_m[wi + 3]);
v = (float4)(matrix_v[wi], matrix_v[wi + 1], matrix_v[wi + 2], matrix_v[wi + 3]);
break;
}
float4 g = matrix_g[i] * inp;
float4 mt = b1 * m + (1 - b1) * g;
float4 vt = sqrt(b2 * v + (1 - b2) * pow(g, 2.0f));
float4 delta = clamp(l * mt / (vt > 0 ? vt : l * 10) - l * WEIGHT_DECAY * weight, -MAX_WEIGHT_DELTA, MAX_WEIGHT_DELTA);
// See UpdateWeightsMomentum's comment just above for why every matrix_w write below is clamped to
// MAX_WEIGHT - this dense Adam kernel is the one that updates the classification output layer.
// No sign-agreement gate on the matrix_w writes (removed 2026-07, all four backends in sync) -
// see AI\NeuronCPU.mqh's CNeuron::updateInputWeights comment for why: it rectified the one-hot
// softmax-CCE gradient stream into a permanent downward ratchet (the all-Neutral collapse).
switch(inputs - j * 4)
{
case 2:
matrix_w[wi + 2] = clamp(matrix_w[wi + 2] + delta.s2, -MAX_WEIGHT, MAX_WEIGHT);
matrix_m[wi + 2] = mt.s2;
matrix_v[wi + 2] = vt.s2;
case 1:
matrix_w[wi + 1] = clamp(matrix_w[wi + 1] + delta.s1, -MAX_WEIGHT, MAX_WEIGHT);
matrix_m[wi + 1] = mt.s1;
matrix_v[wi + 1] = vt.s1;
case 0:
matrix_w[wi] = clamp(matrix_w[wi] + delta.s0, -MAX_WEIGHT, MAX_WEIGHT);
matrix_m[wi] = mt.s0;
matrix_v[wi] = vt.s0;
break;
default:
matrix_w[wi] = clamp(matrix_w[wi] + delta.s0, -MAX_WEIGHT, MAX_WEIGHT);
matrix_m[wi] = mt.s0;
matrix_v[wi] = vt.s0;
matrix_w[wi + 1] = clamp(matrix_w[wi + 1] + delta.s1, -MAX_WEIGHT, MAX_WEIGHT);
matrix_m[wi + 1] = mt.s1;
matrix_v[wi + 1] = vt.s1;
matrix_w[wi + 2] = clamp(matrix_w[wi + 2] + delta.s2, -MAX_WEIGHT, MAX_WEIGHT);
matrix_m[wi + 2] = mt.s2;
matrix_v[wi + 2] = vt.s2;
matrix_w[wi + 3] = clamp(matrix_w[wi + 3] + delta.s3, -MAX_WEIGHT, MAX_WEIGHT);
matrix_m[wi + 3] = mt.s3;
matrix_v[wi + 3] = vt.s3;
break;
}
};
//+------------------------------------------------------------------+
//| Max-pooling layer (CNeuronPoolOCL) - no weights, just a sliding |
//| max over the previous layer's output. Ported from the NeuroNet_DNG|
//| reference's CNeuronProofOCL kernels (these two had no bug). |
//+------------------------------------------------------------------+
__kernel void FeedForwardProof(__global float *matrix_i,
__global float *matrix_o,
int inputs, int window, int step)
{
int i = get_global_id(0);
int pos = i * step;
// Mirrors the bounds guard in this kernel's pure-MQL5 host counterpart
// (CNeuronPoolOCL::feedForwardCPU in AI\NeuronOCLConvPool.mqh) - outputs is sized so pos shouldn't
// reach inputs in practice, but without this guard an out-of-bounds pos would read past matrix_i.
if(pos >= inputs)
return;
float result = matrix_i[pos];
for(int k = 1; k < window; k++)
{
int shift = k + pos;
if(shift >= inputs)
break;
result = max(result, matrix_i[shift]);
}
matrix_o[i] = result;
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
__kernel void CalcInputGradientProof(__global float *matrix_i,
__global float *matrix_g,
__global float *matrix_o,
__global float *matrix_ig,
int outputs, int window, int step)
{
int i = get_global_id(0);
float prev_gradient = 0.0f;
float value = matrix_i[i];
int start = i - window + step;
start = (start - start % step) / step;
int stop = (i - i % step) / step + 1;
for(int out = max(0, start); out < min(outputs, stop); out++)
if(value == matrix_o[out])
prev_gradient += matrix_g[out];
matrix_ig[i] = prev_gradient;
}
//+------------------------------------------------------------------+
//| Convolution layer (CNeuronConvOCL) - ported from the NeuroNet_DNG |
//| reference library and adapted to this project's buffer convention.|
//| window_in+1 weights (incl. bias) are shared across all sliding |
//| positions; window_out lets one position produce several output |
//| channels (filters) at once. |
//+------------------------------------------------------------------+
__kernel void FeedForwardConv(__global float *matrix_w,
__global float *matrix_i,
__global float *matrix_o,
int inputs, int step, int window_in, int window_out,
int activation)
{
int i = get_global_id(0);
int shift_out = window_out * i;
int shift_in = step * i;
for(int out = 0; out < window_out; out++)
{
int shift = (window_in + 1) * out;
int stop = (window_in <= (inputs - shift_in) ? window_in : (inputs - shift_in));
float sum = 0.0f;
for(int k = 0; k < stop; k++)
sum += matrix_i[shift_in + k] * matrix_w[shift + k];
sum += matrix_w[shift + window_in]; // bias
switch(activation)
{
case 0:
sum = tanh(sum);
break;
case 1:
sum = 1 / (1 + exp(-clamp(sum, -50.0f, 50.0f)));
break;
case 2: // PReLU, param=0.01
if(sum < 0)
sum *= 0.01f;
break;
}
matrix_o[out + shift_out] = sum;
}
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
__kernel void CalcHiddenGradientConv(__global float *matrix_w,
__global float *matrix_g,
__global float *matrix_o,
__global float *matrix_ig,
int outputs, int step, int window_in, int window_out,
int activation)
{
int i = get_global_id(0);
float out = matrix_o[i];
int start = i - window_in + step;
start = max((start - start % step) / step, 0);
int stop = (i - i % step) / step + 1;
if(stop > (outputs / window_out))
stop = outputs / window_out;
float sum = 0.0f;
for(int h = 0; h < window_out; h++)
{
for(int k = start; k < stop; k++)
{
int shift_w = (stop - k - 1) * step + i % step + h * (window_in + 1);
int shift_g = k * window_out + h;
if(shift_g >= outputs || shift_w >= (window_in + 1) * window_out)
break;
sum += matrix_g[shift_g] * matrix_w[shift_w];
}
}
switch(activation)
{
case 0:
sum = clamp(sum + out, -1.0f, 1.0f) - out;
sum = sum * fmax(MIN_ACTIVATION_DERIVATIVE, 1 - pow(out, 2.0f));
break;
case 1:
sum = clamp(sum + out, 0.0f, 1.0f) - out;
sum = sum * fmax(MIN_ACTIVATION_DERIVATIVE, out * (1 - out));
break;
case 2: // PReLU, param=0.01
if(out < 0)
sum *= 0.01f;
break;
}
matrix_ig[i] = sum;
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
__kernel void UpdateWeightsConvMomentum(__global float *matrix_w,
__global float *matrix_g,
__global float *matrix_i,
__global float *matrix_dw,
int inputs, float learning_rates, float momentum,
int window_in, int window_out, int step, int optimizer)
{
int i = get_global_id(0);
int shift = i % (window_in + 1);
int shift_out = (i - shift) / (window_in + 1);
int total = (inputs - window_in) % step;
total = (inputs - window_in - total) / step + (total > 0 ? 1 : 0);
float grad = 0.0f;
for(int t = 0; t < total; t++)
{
if(shift != window_in && (shift + t * step) >= inputs)
break;
grad += matrix_g[t * window_out + shift_out] * (shift == window_in ? 1 : matrix_i[shift + t * step]);
}
// `optimizer` retained unused - see UpdateWeightsMomentum's comment above.
float delta = learning_rates * grad + momentum * matrix_dw[i];
matrix_dw[i] = delta;
matrix_w[i] = clamp(matrix_w[i] + delta, -MAX_WEIGHT, MAX_WEIGHT);
};
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
__kernel void UpdateWeightsConvAdam(__global float *matrix_w,
__global const float *matrix_g,
__global const float *matrix_i,
__global float *matrix_m,
__global float *matrix_v,
const int inputs, const float l, const float b1, const float b2,
int window_in, int window_out, int step)
{
int i = get_global_id(0); // weight offset within a single output channel's (window_in+1) block
if(i > window_in)
return;
int total = (inputs - (window_in - step)) % step;
total = (inputs - (window_in - step) - total) / step + (total > 0 ? 1 : 0);
for(int out = 0; out < window_out; out++)
{
int shift_w = i + out * (window_in + 1);
float grad = 0.0f;
for(int t = 0; t < total; t++)
{
if(i != window_in && (i + t * step) >= inputs)
break;
grad += matrix_g[t * window_out + out] * (i == window_in ? 1 : matrix_i[i + t * step]);
}
float mt = b1 * matrix_m[shift_w] + (1 - b1) * grad;
float vt = sqrt(b2 * matrix_v[shift_w] + (1 - b2) * pow(grad, 2.0f));
float delta = clamp(l * mt / (vt > 0 ? vt : l * 10) - l * WEIGHT_DECAY * matrix_w[shift_w], -MAX_WEIGHT_DELTA, MAX_WEIGHT_DELTA);
// Sign-agreement gate removed - see UpdateWeightsAdam's comment above for why.
matrix_w[shift_w] = clamp(matrix_w[shift_w] + delta, -MAX_WEIGHT, MAX_WEIGHT);
matrix_m[shift_w] = mt;
matrix_v[shift_w] = vt;
}
};
//+------------------------------------------------------------------+
//| LSTM (CNeuronLSTMOCL) - derived from scratch from the standard |
//| LSTM equations (NOT ported from any reference implementation). |
//| Single-timestep-truncated BPTT: gradient does not flow back into |
//| h_prev/c_prev from a previous step - only this step's gates and |
//| weights receive gradient. Adam-only optimizer. |
//| Weight layout: 4 gates (order: forget, input, output, candidate), |
//| each a block of H*(I+H+1) floats - per hidden unit: H weights on |
//| h_prev, I weights on the current input, 1 bias. |
//+------------------------------------------------------------------+
__kernel void LSTM_Gates(__global float *matrix_w,
__global float *hidden_prev,
__global float *inputs,
__global float *concatenated,
int hidden_size, int input_size)
{
int id = get_global_id(0);
int gate = get_global_id(1);
int per_gate = hidden_size * (hidden_size + input_size + 1);
int shift = gate * per_gate + id * (hidden_size + input_size + 1);
float sum = 0.0f;
for(int k = 0; k < hidden_size; k++)
sum += hidden_prev[k] * matrix_w[shift + k];
for(int k = 0; k < input_size; k++)
sum += inputs[k] * matrix_w[shift + hidden_size + k];
sum += matrix_w[shift + hidden_size + input_size];
float val = (gate < 3) ? (1.0f / (1.0f + exp(-sum))) : tanh(sum);
concatenated[gate * hidden_size + id] = val;
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
__kernel void LSTM_State(__global float *concatenated,
__global float *memory,
__global float *hidden_prev,
__global float *hidden_cache,
__global float *output,
int hidden_size)
{
int id = get_global_id(0);
float f = concatenated[id];
float ii = concatenated[hidden_size + id];
float o = concatenated[2 * hidden_size + id];
float g = concatenated[3 * hidden_size + id];
float c_prev = memory[id];
memory[hidden_size + id] = c_prev;
float c_t = f * c_prev + ii * g;
memory[id] = c_t;
hidden_cache[id] = hidden_prev[id];
output[id] = o * tanh(c_t);
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
__kernel void LSTM_GateGradient(__global float *gradient,
__global float *memory,
__global float *concatenated,
__global float *concatenated_gradient,
int hidden_size)
{
int id = get_global_id(0);
float c_t = memory[id];
float c_prev = memory[hidden_size + id];
float f = concatenated[id];
float ii = concatenated[hidden_size + id];
float o = concatenated[2 * hidden_size + id];
float g = concatenated[3 * hidden_size + id];
float t = tanh(c_t);
float dh = gradient[id];
// Derivative floor on every gate below - see MIN_ACTIVATION_DERIVATIVE's declaration comment.
// LSTM gates are sigmoid/tanh by design and routinely sit near their saturated extremes (a fully
// "open" or "closed" gate), which is exactly where an unfloored derivative goes to zero and
// blocks gradient from flowing through the recurrent path at all.
float dc = dh * o * fmax(MIN_ACTIVATION_DERIVATIVE, 1 - t * t);
concatenated_gradient[2 * hidden_size + id] = dh * t * fmax(MIN_ACTIVATION_DERIVATIVE, o * (1 - o)); // output gate
concatenated_gradient[id] = dc * c_prev * fmax(MIN_ACTIVATION_DERIVATIVE, f * (1 - f)); // forget gate
concatenated_gradient[hidden_size + id] = dc * g * fmax(MIN_ACTIVATION_DERIVATIVE, ii * (1 - ii)); // input gate
concatenated_gradient[3 * hidden_size + id] = dc * ii * fmax(MIN_ACTIVATION_DERIVATIVE, 1 - g * g); // candidate
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
__kernel void LSTM_WeightsGradient(__global float *concatenated_gradient,
__global float *hidden_cache,
__global float *inputs,
__global float *weights_gradient,
int hidden_size, int input_size)
{
int wi = get_global_id(0);
int per_gate = hidden_size * (hidden_size + input_size + 1);
int gate = wi / per_gate;
int rem = wi % per_gate;
int id = rem / (hidden_size + input_size + 1);
int k = rem % (hidden_size + input_size + 1);
float inp = (k < hidden_size ? hidden_cache[k] : (k < hidden_size + input_size ? inputs[k - hidden_size] : 1.0f));
weights_gradient[wi] = concatenated_gradient[gate * hidden_size + id] * inp;
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
__kernel void LSTM_InputsGradient(__global float *concatenated_gradient,
__global float *matrix_w,
__global float *inputs_gradient,
int hidden_size, int input_size)
{
int j = get_global_id(0);
int per_gate = hidden_size * (hidden_size + input_size + 1);
float sum = 0.0f;
for(int gate = 0; gate < 4; gate++)
for(int id = 0; id < hidden_size; id++)
sum += concatenated_gradient[gate * hidden_size + id] * matrix_w[gate * per_gate + id * (hidden_size + input_size + 1) + hidden_size + j];
inputs_gradient[j] = sum;
}
//+------------------------------------------------------------------+
//| |
//+------------------------------------------------------------------+
__kernel void LSTM_UpdateWeightsAdam(__global float *matrix_w,
__global float *weights_gradient,
__global float *matrix_m,
__global float *matrix_v,
float l, float b1, float b2)
{
int wi = get_global_id(0);
float g = weights_gradient[wi];
float mt = b1 * matrix_m[wi] + (1 - b1) * g;
float vt = sqrt(b2 * matrix_v[wi] + (1 - b2) * g * g);
float delta = clamp(l * mt / (vt > 0 ? vt : l * 10) - l * WEIGHT_DECAY * matrix_w[wi], -MAX_WEIGHT_DELTA, MAX_WEIGHT_DELTA);
// Sign-agreement gate removed - see UpdateWeightsAdam's comment above for why.
matrix_w[wi] = clamp(matrix_w[wi] + delta, -MAX_WEIGHT, MAX_WEIGHT);
matrix_m[wi] = mt;
matrix_v[wi] = vt;
};
//+------------------------------------------------------------------+
//| SGD+momentum counterpart to LSTM_UpdateWeightsAdam above - same |
//| flat (already-elementwise) weights_gradient/matrix_dw layout, just|
//| the classic heavy-ball update instead of Adam's per-parameter step|
//+------------------------------------------------------------------+
__kernel void LSTM_UpdateWeightsMomentum(__global float *matrix_w,
__global float *weights_gradient,
__global float *matrix_dw,
float learning_rates, float momentum, int optimizer)
{
int wi = get_global_id(0);
// `optimizer` retained unused - see UpdateWeightsMomentum's comment above.
float delta = learning_rates * weights_gradient[wi] + momentum * matrix_dw[wi];
matrix_dw[wi] = delta;
matrix_w[wi] = clamp(matrix_w[wi] + delta, -MAX_WEIGHT, MAX_WEIGHT);
};