Summary

While testing hls4ml 1.3.0, i noticed inconsistent classification results when converting models that use a softmax output layer.
converting a model with a softmax output activation to HLS C-simulation via hls_model.compile() + hls_model.predict() causes a large fraction of argmax classification decisions to differ from the float32 Keras reference, despite producing correctly normalized probability outputs (Σp = 1 in all runs).

In my experiments, the number of mismatched predictions ranged from 0 to 200 out of 200 samples, depending on the model initialization.
To isolate the source of the error, i also converted the same model without the final softmax layer and applied softmax in NumPy after HLS inference. Doing this reduced the mismatches from 494 total flips to 18 across the same eight test runs, while using identical weights and inputs ....this suggests the issue is a LUT approximation error.

Relation to #1443: Issue #1443 reports that the latency softmax
implementation fails tests and produces ~19% accuracy,this issue
provides a controlled reproduction using the default softmax strategy
(not latency), isolates the LUT as the sole cause via a logit only
control experiment, and establishes that the failure is weight distribution
dependent ,absent for some random seeds, severe (up to 100% argmax flips)
for others, the normalization is always correct (Σpk = 1.0), making
this a silent failure with no error raised.

Environment

hls4ml 1.3.0
TensorFlow 2.19.0
tf_keras 2.19.0
Python 3.12
Platform Google Colab (Linux x86_64)
FPGA target - xc7z020clg400-1 (Vivado backend)
Precision - ap_fixed<24,12>, ReuseFactor=1, Latency strategy

Minimal Reproducible Example

import numpy as np
import tf_keras as keras
import hls4ml, shutil

import keras as _k
_k.__version__ = "2.15.0"

np.random.seed(0)
N, T = 200, 140
X = np.random.randn(N, T, 1).astype(np.float32)

def numpy_softmax(logits):
    e = np.exp(logits - logits.max(axis=1, keepdims=True))
    return e / e.sum(axis=1, keepdims=True)

def build_and_test(weight_seed, precision="ap_fixed<24,12>"):
    # full model with softmax
    keras.utils.set_random_seed(weight_seed)
    inp    = keras.Input(shape=(T, 1))
    x      = keras.layers.GlobalAveragePooling1D()(inp)
    x      = keras.layers.Dense(8, activation="relu")(x)
    out    = keras.layers.Dense(2, activation="softmax")(x)
    full_m = keras.Model(inp, out)

    # logit model with identical weights, no softmax
    keras.utils.set_random_seed(weight_seed)
    inp2    = keras.Input(shape=(T, 1))
    x2      = keras.layers.GlobalAveragePooling1D()(inp2)
    x2      = keras.layers.Dense(8, activation="relu")(x2)
    out2    = keras.layers.Dense(2)(x2)
    logit_m = keras.Model(inp2, out2)
    for i in range(len(logit_m.layers)):
        try: logit_m.layers[i].set_weights(full_m.layers[i].get_weights())
        except: pass

    keras_pred = np.argmax(full_m.predict(X, verbose=0), axis=1)

    def hls_predict(model, out_dir, is_logit=False):
        shutil.rmtree(out_dir, ignore_errors=True)
        cfg = hls4ml.utils.config_from_keras_model(
            model, granularity="name", default_precision=precision)
        for l in cfg["LayerName"]:
            cfg["LayerName"][l]["Strategy"]    = "Latency"
            cfg["LayerName"][l]["ReuseFactor"] = 1
        hls_m = hls4ml.converters.convert_from_keras_model(
            model, hls_config=cfg,
            output_dir=out_dir, part="xc7z020clg400-1")
        hls_m.compile()
        raw = hls_m.predict(np.ascontiguousarray(X, dtype=np.float32))
        return numpy_softmax(raw) if is_logit else raw

    hls_full  = hls_predict(full_m,  f"hls_full_{weight_seed}")
    hls_logit = hls_predict(logit_m, f"hls_logit_{weight_seed}", is_logit=True)

    return {
        "flips_full":  int((np.argmax(hls_full,  axis=1) != keras_pred).sum()),
        "flips_logit": int((np.argmax(hls_logit, axis=1) != keras_pred).sum()),
        "sum_err":     float(np.abs(hls_full.sum(axis=1) - 1.0).max()),
    }

for seed in [0, 1, 7, 42, 99, 123, 256, 999]:
    r = build_and_test(seed)
    print(f"Seed {seed:>4}: HLS softmax {r['flips_full']:>3}/200  "
          f"Logit-only {r['flips_logit']:>3}/200  "
          f"max|Σpk-1|={r['sum_err']:.6f}")

Observed Output

Seed    0: HLS softmax  94/200   Logit-only   0/200   max|Σpk-1|=0.000000
Seed    1: HLS softmax 200/200   Logit-only  10/200   max|Σpk-1|=0.000000
Seed    7: HLS softmax   0/200   Logit-only   3/200   max|Σpk-1|=0.000000
Seed   42: HLS softmax   0/200   Logit-only   1/200   max|Σpk-1|=0.000000
Seed   99: HLS softmax  94/200   Logit-only   0/200   max|Σpk-1|=0.000000
Seed  123: HLS softmax   0/200   Logit-only   0/200   max|Σpk-1|=0.000000
Seed  256: HLS softmax   0/200   Logit-only   0/200   max|Σpk-1|=0.000000
Seed  999: HLS softmax 106/200   Logit-only   4/200   max|Σpk-1|=0.000000

Possible Cause

From tracing through the implementation, my current understanding is that the issue is likely related to the fixed-point exponential lookup table used by the HLS softmax implementation (nnet_activation.h),.
When two logits are very close, small approximation errors in the exponential values appear sufficient to change which probability becomes slightly larger after normalization. The probabilities still sum to one, but the predicted class changes.

Workaround

For now, I was able to avoid the issue by converting the model without the final softmax layer and applying softmax in Python after HLS inference.

def numpy_softmax(logits):
    e = np.exp(logits - logits.max(axis=1, keepdims=True))
    return e / e.sum(axis=1, keepdims=True)

hw_prob = numpy_softmax(hls_logit_model.predict(X))

With the same weights and inputs, this reduced the worst-case mismatch from 200/200 predictions to 10/200.

Why this matters

This affects any hls4ml user who:

1. Has a softmax output activation in their model
2. Uses hls_model.compile() + hls_model.predict() for C-simulation
3. Has a model that produces near-boundary predictions on some inputs

The failure is entirely silent , valid looking normalized probabilities, no error, no warning.

Possible Improvements

It might be useful to:

1. Investigate whether the exponential LUT precision in nnet_activation.h can be improved.
2. Add a regression test that checks argmax agreement, not only probability normalization.
3. Document this behavior or suggest converting the logit model when exact agreement with the original model is important.

Reference

I ran into this issue while deploying an SE-attention network for ECG anomaly
detection on a Xilinx PYNQ-Z2 board using hls4ml.

The full reproduction code, along with the validation experiments discussed
above, is available here:

https://github.com/niteshg97/ecg-fpga-se-networks