+ New

utensil

Public
120b3201f366400dfe818461ee828002bf8ef9ca
diff --git a/.gitignore b/.gitignore
index 672c3aa..167863f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,8 +4,10 @@
 *.out
 test
 mnist/
+speech_commands/
 examples/mnist_mlp
 examples/mnist_cnn
 examples/mobilenetv1_cifar10
 examples/resnet8_cifar10
+examples/lstm_kws
 
diff --git a/Makefile b/Makefile
index c82f5e1..98106e3 100644
--- a/Makefile
+++ b/Makefile
@@ -31,9 +31,20 @@ cifar/data_batch_1.bin:
 
 cifar-data: cifar/data_batch_1.bin
 
+SC_URL := http://download.tensorflow.org/data/speech_commands_v0.02.tar.gz
+
+speech_commands/testing_list.txt:
+	mkdir -p speech_commands
+	curl -sL $(SC_URL) -o speech_commands/speech_commands_v0.02.tar.gz
+	tar xzf speech_commands/speech_commands_v0.02.tar.gz -C speech_commands
+	rm -f speech_commands/speech_commands_v0.02.tar.gz
+
+speech-commands-data: speech_commands/testing_list.txt
+
 EXAMPLES := $(patsubst %.c,%,$(wildcard examples/*.c))
 MNIST_EXAMPLES := examples/mnist_mlp examples/mnist_cnn
 CIFAR_EXAMPLES := examples/resnet8_cifar10 examples/mobilenetv1_cifar10
+LSTM_EXAMPLES := examples/lstm_kws
 
 examples: $(EXAMPLES)
 
@@ -42,5 +53,6 @@ examples/%: examples/%.c utensil.h
 
 $(MNIST_EXAMPLES): examples/mnist_loader.h mnist-data
 $(CIFAR_EXAMPLES): examples/cifar10_loader.h cifar-data
+$(LSTM_EXAMPLES): examples/speech_commands_loader.h speech-commands-data
 
-.PHONY: test examples mnist-data cifar-data
+.PHONY: test examples mnist-data cifar-data speech-commands-data
diff --git a/test.c b/test.c
index 83af798..47accfb 100644
--- a/test.c
+++ b/test.c
@@ -451,6 +451,148 @@ static void test_conv1d_backward(ut_dev dev) {
   ut_conv1d_free(&l);
 }
 
+static void test_lstm_step(ut_dev dev) {
+  // in=1, hidden=1, B=1. Values chosen so every gate sits away from 0/1
+  // saturation; expected numbers cross-checked against a reference
+  // double-precision implementation of the same formulas.
+  ut_lstm l = ut_lstm_alloc(1, 1, dev);
+  ut_free(l.W_ih);
+  l.W_ih = ut_from_data(2, (int[]){1, 4}, (float[]){1.f, 0.5f, -1.f, 0.3f}, dev);
+  ut_free(l.W_hh);
+  l.W_hh = ut_from_data(2, (int[]){1, 4}, (float[]){0.2f, -0.3f, 0.4f, -0.1f}, dev);
+  ut_free(l.bias);
+  l.bias = ut_from_data(1, (int[]){4}, (float[]){0.1f, 0.2f, -0.2f, 0.05f}, dev);
+
+  ut_tensor* x = ut_from_data(2, (int[]){1, 1}, (float[]){1.f}, dev);
+  ut_tensor* h_prev = ut_from_data(2, (int[]){1, 1}, (float[]){0.5f}, dev);
+  ut_tensor* c_prev = ut_from_data(2, (int[]){1, 1}, (float[]){0.2f}, dev);
+
+  ut_lstm_cache cache;
+  ut_tensor* c_out;
+  ut_tensor* h = ut_lstm_step(&l, x, h_prev, c_prev, &c_out, &cache);
+  ut_sync_cpu(h);
+  ut_sync_cpu(c_out);
+  assert_eq(h->data[0], -0.24634508f, 1e-4f);
+  assert_eq(c_out->data[0], -0.45847687f, 1e-4f);
+
+  ut_tensor* dh = ut_from_data(2, (int[]){1, 1}, (float[]){1.f}, dev);
+  ut_tensor* dc = ut_alloc(2, (int[]){1, 1}, dev);  // zero: c only feeds the next step
+  ut_tensor* dW_ih = ut_alloc(2, (int[]){1, 4}, dev);
+  ut_tensor* dW_hh = ut_alloc(2, (int[]){1, 4}, dev);
+  ut_tensor* db = ut_alloc(1, (int[]){4}, dev);
+  ut_tensor *dx, *dh_prev, *dc_prev;
+  ut_lstm_backward(&l, &cache, dh, dc, dW_ih, dW_hh, db, &dx, &dh_prev, &dc_prev);
+  ut_sync_cpu(dW_ih);
+  ut_sync_cpu(dW_hh);
+  ut_sync_cpu(db);
+  ut_sync_cpu(dx);
+  ut_sync_cpu(dh_prev);
+  ut_sync_cpu(dc_prev);
+
+  assert_data(dW_ih, ((float[]){-0.06351452f, 0.02175301f, 0.15131002f, -0.10483399f}), 1e-4f);
+  assert_data(dW_hh, ((float[]){-0.03175726f, 0.01087650f, 0.07565501f, -0.05241700f}), 1e-4f);
+  assert_data(db, ((float[]){-0.06351452f, 0.02175301f, 0.15131002f, -0.10483399f}), 1e-4f);
+  assert_eq(dx->data[0], -0.23539823f, 1e-4f);
+  assert_eq(dh_prev->data[0], 0.05177860f, 1e-4f);
+  assert_eq(dc_prev->data[0], 0.29728238f, 1e-4f);
+
+  ut_lstm_cache_free(&cache);
+  ut_free_all(x, h_prev, c_prev, h, c_out, dh, dc, dW_ih, dW_hh, db, dx, dh_prev, dc_prev);
+  ut_lstm_free(&l);
+}
+
+static float _lstm_seq_loss(ut_lstm* l, ut_tensor* x_seq) {
+  ut_lstm_seq_cache cache;
+  ut_tensor *h_n, *c_n;
+  ut_tensor* h_seq = ut_lstm_forward_seq(l, x_seq, NULL, NULL, &cache, &h_n, &c_n);
+  ut_sync_cpu(h_n);
+  float loss = 0.f;
+  for (int k = 0; k < h_n->shape.nelem; k++) loss += h_n->data[k];
+  ut_lstm_seq_cache_free(&cache);
+  ut_free_all(h_seq, h_n, c_n);
+  return loss;
+}
+
+// Sequence-level BPTT is exercised via a numerical gradient check (loss =
+// sum of the final hidden state) rather than hand algebra, since a 3-step
+// unroll is impractical to derive by hand.
+static void test_lstm_seq_gradcheck(void) {
+  ut_dev dev = UT_CPU;
+  int T = 3, B = 1, IN = 2, H = 2;
+  ut_lstm l = ut_lstm_alloc(IN, H, dev);
+  ut_free(l.W_ih);
+  l.W_ih = ut_from_data(
+      2, (int[]){IN, 4 * H},
+      (float[]){0.3f, -0.2f, 0.1f, 0.4f, -0.1f, 0.5f, 0.2f, -0.3f, 0.15f, 0.1f, -0.2f, 0.05f,
+                0.25f, -0.15f, 0.1f, 0.2f},
+      dev);
+  ut_free(l.W_hh);
+  l.W_hh = ut_from_data(
+      2, (int[]){H, 4 * H},
+      (float[]){0.15f, -0.25f, 0.05f, 0.2f, -0.1f, 0.3f, -0.2f, 0.1f, 0.05f, 0.1f, -0.15f, 0.2f,
+                -0.05f, 0.25f, 0.1f, -0.1f},
+      dev);
+  ut_free(l.bias);
+  l.bias = ut_from_data(1, (int[]){4 * H}, (float[]){0.1f, 0.2f, -0.1f, 0.05f, 0.f, -0.05f, 0.1f, 0.f},
+                        dev);
+
+  float xdata[3 * 1 * 2] = {0.5f, -0.3f, 0.2f, 0.1f, -0.4f, 0.6f};
+  ut_tensor* x_seq = ut_from_data(3, (int[]){T, B, IN}, xdata, dev);
+
+  ut_lstm_seq_cache cache;
+  ut_tensor *h_n, *c_n;
+  ut_tensor* h_seq = ut_lstm_forward_seq(&l, x_seq, NULL, NULL, &cache, &h_n, &c_n);
+
+  ut_tensor* grad_h_seq = ut_alloc(3, (int[]){T, B, H}, dev);  // zero except last step
+  for (int k = 0; k < B * H; k++) grad_h_seq->data[(T - 1) * B * H + k] = 1.f;
+
+  ut_tensor* dW_ih = ut_alloc(2, (int[]){IN, 4 * H}, dev);
+  ut_tensor* dW_hh = ut_alloc(2, (int[]){H, 4 * H}, dev);
+  ut_tensor* db = ut_alloc(1, (int[]){4 * H}, dev);
+  ut_tensor* dx_seq = ut_lstm_backward_seq(&l, &cache, grad_h_seq, dW_ih, dW_hh, db);
+  ut_sync_cpu(dW_ih);
+  ut_sync_cpu(dW_hh);
+  ut_sync_cpu(dx_seq);
+
+  ut_lstm_seq_cache_free(&cache);
+  ut_free_all(h_seq, h_n, c_n, grad_h_seq);
+
+  float eps = 1e-3f;
+  {
+    int idx = 3;
+    float orig = l.W_ih->data[idx];
+    l.W_ih->data[idx] = orig + eps;
+    float lp = _lstm_seq_loss(&l, x_seq);
+    l.W_ih->data[idx] = orig - eps;
+    float lm = _lstm_seq_loss(&l, x_seq);
+    l.W_ih->data[idx] = orig;
+    assert_eq(dW_ih->data[idx], (lp - lm) / (2.f * eps), 2e-2f);
+  }
+  {
+    int idx = 5;
+    float orig = l.W_hh->data[idx];
+    l.W_hh->data[idx] = orig + eps;
+    float lp = _lstm_seq_loss(&l, x_seq);
+    l.W_hh->data[idx] = orig - eps;
+    float lm = _lstm_seq_loss(&l, x_seq);
+    l.W_hh->data[idx] = orig;
+    assert_eq(dW_hh->data[idx], (lp - lm) / (2.f * eps), 2e-2f);
+  }
+  {
+    int idx = 2;  // t=1, b=0, in=0
+    float orig = x_seq->data[idx];
+    x_seq->data[idx] = orig + eps;
+    float lp = _lstm_seq_loss(&l, x_seq);
+    x_seq->data[idx] = orig - eps;
+    float lm = _lstm_seq_loss(&l, x_seq);
+    x_seq->data[idx] = orig;
+    assert_eq(dx_seq->data[idx], (lp - lm) / (2.f * eps), 2e-2f);
+  }
+
+  ut_free_all(x_seq, dW_ih, dW_hh, db, dx_seq);
+  ut_lstm_free(&l);
+}
+
 static void test_conv2d_forward(ut_dev dev) {
   // N=2,in_c=1,out_c=2,kh=kw=2,s=1,p=0. Distinct per-batch/per-channel bias
   // catches any N<->C mixup in the forward transpose or bias broadcast.
@@ -733,6 +875,12 @@ int main() {
   test_conv1d_backward(UT_METAL);
 #endif
 
+  test_lstm_step(UT_CPU);
+#ifdef __APPLE__
+  test_lstm_step(UT_METAL);
+#endif
+  test_lstm_seq_gradcheck();
+
   test_conv2d_forward(UT_CPU);
 #ifdef __APPLE__
   test_conv2d_forward(UT_METAL);
diff --git a/utensil.h b/utensil.h
index 4657e0a..cc32433 100644
--- a/utensil.h
+++ b/utensil.h
@@ -59,6 +59,25 @@ typedef struct ut_conv1d_cache {
   ut_tensor* col;
 } ut_conv1d_cache;
 
+typedef struct ut_lstm {
+  ut_tensor* W_ih;  // [in, 4*hidden] gate order: i, f, g, o
+  ut_tensor* W_hh;  // [hidden, 4*hidden]
+  ut_tensor* bias;  // [4*hidden]
+  int in, hidden;
+} ut_lstm;
+
+// One timestep's worth of state, retained for that step's backward pass.
+// Gate math always runs on CPU (see ut_lstm_step), so these are CPU tensors.
+typedef struct ut_lstm_cache {
+  ut_tensor *x, *h_prev, *c_prev;  // retained inputs to this step
+  ut_tensor *i, *f, *g, *o, *c;    // gate activations and new cell state
+} ut_lstm_cache;
+
+typedef struct ut_lstm_seq_cache {
+  ut_lstm_cache* steps;  // [t]
+  int t;
+} ut_lstm_seq_cache;
+
 typedef struct ut_conv2d {
   ut_tensor* weight;  // [out_c, in_c, kh, kw]
   ut_tensor* bias;    // [out_c]
@@ -1754,6 +1773,250 @@ void ut_conv1d_free(ut_conv1d* l) {
   ut_free_all(l->weight, l->bias);
   l->weight = l->bias = NULL;
 }
+
+// =========================================================
+// LSTM
+// =========================================================
+// Vanilla single-layer LSTM. The two big matmuls per step (x@W_ih, h@W_hh)
+// go through ut_matmul, so they run on whichever device the weights/inputs
+// are on; the gating nonlinearities and state update are a plain CPU loop,
+// the same way conv1d/conv2d finish their bias-add and layernorm/batchnorm2d
+// compute their statistics. Stacking layers or running bidirectionally is
+// left to the caller — wire two ut_lstm instances together.
+
+ut_lstm ut_lstm_alloc(int in, int hidden, ut_dev dev) {
+  ut_lstm l = {.in = in, .hidden = hidden};
+  float std = sqrtf(1.f / (float)hidden);
+  l.W_ih = ut_randn(2, (int[]){in, 4 * hidden}, 0.f, std, dev);
+  l.W_hh = ut_randn(2, (int[]){hidden, 4 * hidden}, 0.f, std, dev);
+  l.bias = ut_alloc(1, (int[]){4 * hidden}, dev);
+  // forget-gate bias = 1 (Jozefowicz et al. 2015) so early training doesn't
+  // forget everything by default
+  ut_sync_cpu(l.bias);
+  for (int j = 0; j < hidden; j++) l.bias->data[hidden + j] = 1.f;
+  l.bias->dirty_gpu = true;
+  return l;
+}
+
+// x:[B,in], h_prev/c_prev:[B,hidden] -> returns h:[B,hidden], writes new cell
+// state into *c_out:[B,hidden] (required; caller owns the returned tensor and
+// typically feeds it back in as c_prev on the next step).
+ut_tensor* ut_lstm_step(ut_lstm* l, ut_tensor* x, ut_tensor* h_prev, ut_tensor* c_prev,
+                        ut_tensor** c_out, ut_lstm_cache* cache) {
+  int B = x->shape.shape[0], H = l->hidden;
+
+  ut_tensor* zx = ut_matmul(x, l->W_ih);       // [B,4H]
+  ut_tensor* zh = ut_matmul(h_prev, l->W_hh);  // [B,4H]
+  ut_tensor* z = ut_add(zx, zh);
+  ut_free_all(zx, zh);
+  ut_sync_cpu(z);
+  ut_sync_cpu(l->bias);
+  ut_sync_cpu(c_prev);
+
+  ut_tensor* gi = ut_alloc(2, (int[]){B, H}, UT_CPU);
+  ut_tensor* gf = ut_alloc(2, (int[]){B, H}, UT_CPU);
+  ut_tensor* gg = ut_alloc(2, (int[]){B, H}, UT_CPU);
+  ut_tensor* go = ut_alloc(2, (int[]){B, H}, UT_CPU);
+  ut_tensor* c = ut_alloc(2, (int[]){B, H}, UT_CPU);
+  ut_tensor* h = ut_alloc(2, (int[]){B, H}, UT_CPU);
+
+  for (int b = 0; b < B; b++) {
+    const float* zb = z->data + b * 4 * H;
+    const float* bb = l->bias->data;
+    const float* cp = c_prev->data + b * H;
+    for (int j = 0; j < H; j++) {
+      float ig = 1.f / (1.f + expf(-(zb[j] + bb[j])));
+      float fg = 1.f / (1.f + expf(-(zb[H + j] + bb[H + j])));
+      float cg = tanhf(zb[2 * H + j] + bb[2 * H + j]);
+      float og = 1.f / (1.f + expf(-(zb[3 * H + j] + bb[3 * H + j])));
+      float cc = fg * cp[j] + ig * cg;
+      gi->data[b * H + j] = ig;
+      gf->data[b * H + j] = fg;
+      gg->data[b * H + j] = cg;
+      go->data[b * H + j] = og;
+      c->data[b * H + j] = cc;
+      h->data[b * H + j] = og * tanhf(cc);
+    }
+  }
+  ut_free(z);
+
+  if (cache) {
+    cache->x = ut_retain(x);
+    cache->h_prev = ut_retain(h_prev);
+    cache->c_prev = ut_retain(c_prev);
+    cache->i = gi;
+    cache->f = gf;
+    cache->g = gg;
+    cache->o = go;
+    cache->c = ut_retain(c);
+  } else {
+    ut_free_all(gi, gf, gg, go);
+  }
+  *c_out = c;
+  return h;
+}
+
+// dh: gradient w.r.t. this step's h (recurrent + any direct consumer).
+// dc: gradient w.r.t. this step's c flowing in from the next step (pass a
+// zeroed [B,hidden] tensor at the last step of a sequence).
+// Accumulates into dW_ih/dW_hh/db (caller-owned accumulators, like every
+// other _backward in this file); writes fresh *dx/*dh_prev/*dc_prev.
+void ut_lstm_backward(ut_lstm* l, ut_lstm_cache* c, ut_tensor* dh, ut_tensor* dc, ut_tensor* dW_ih,
+                      ut_tensor* dW_hh, ut_tensor* db, ut_tensor** dx, ut_tensor** dh_prev,
+                      ut_tensor** dc_prev) {
+  int B = c->x->shape.shape[0], H = l->hidden;
+  ut_sync_cpu(dh);
+  ut_sync_cpu(dc);
+  ut_sync_cpu(c->i);
+  ut_sync_cpu(c->f);
+  ut_sync_cpu(c->g);
+  ut_sync_cpu(c->o);
+  ut_sync_cpu(c->c);
+  ut_sync_cpu(c->c_prev);
+
+  ut_tensor* dz = ut_alloc(2, (int[]){B, 4 * H}, UT_CPU);
+  ut_tensor* dcp = ut_alloc(2, (int[]){B, H}, UT_CPU);
+
+  for (int b = 0; b < B; b++) {
+    const float* dhb = dh->data + b * H;
+    const float* dcb = dc->data + b * H;
+    const float* ib = c->i->data + b * H;
+    const float* fb = c->f->data + b * H;
+    const float* gb = c->g->data + b * H;
+    const float* ob = c->o->data + b * H;
+    const float* cb = c->c->data + b * H;
+    const float* cpb = c->c_prev->data + b * H;
+    float* dzb = dz->data + b * 4 * H;
+    float* dcpb = dcp->data + b * H;
+    for (int j = 0; j < H; j++) {
+      float tc = tanhf(cb[j]);
+      float dov = dhb[j] * tc;
+      float dctot = dcb[j] + dhb[j] * ob[j] * (1.f - tc * tc);
+      float dfv = dctot * cpb[j];
+      float div = dctot * gb[j];
+      float dgv = dctot * ib[j];
+      dcpb[j] = dctot * fb[j];
+      dzb[j] = div * ib[j] * (1.f - ib[j]);
+      dzb[H + j] = dfv * fb[j] * (1.f - fb[j]);
+      dzb[2 * H + j] = dgv * (1.f - gb[j] * gb[j]);
+      dzb[3 * H + j] = dov * ob[j] * (1.f - ob[j]);
+    }
+  }
+
+  ut_tensor* dWx = ut_matmul_t(c->x, dz, true, false);
+  ut_sync_cpu(dWx);
+  ut_sync_cpu(dW_ih);
+  for (int k = 0; k < dW_ih->shape.nelem; k++) dW_ih->data[k] += dWx->data[k];
+  ut_free(dWx);
+
+  ut_tensor* dWh = ut_matmul_t(c->h_prev, dz, true, false);
+  ut_sync_cpu(dWh);
+  ut_sync_cpu(dW_hh);
+  for (int k = 0; k < dW_hh->shape.nelem; k++) dW_hh->data[k] += dWh->data[k];
+  ut_free(dWh);
+
+  if (db) {
+    ut_sync_cpu(db);
+    for (int b = 0; b < B; b++)
+      for (int k = 0; k < 4 * H; k++) db->data[k] += dz->data[b * 4 * H + k];
+  }
+
+  *dx = ut_matmul_t(dz, l->W_ih, false, true);
+  *dh_prev = ut_matmul_t(dz, l->W_hh, false, true);
+  *dc_prev = dcp;
+
+  ut_free(dz);
+}
+
+void ut_lstm_cache_free(ut_lstm_cache* c) {
+  ut_free_all(c->x, c->h_prev, c->c_prev, c->i, c->f, c->g, c->o, c->c);
+  memset(c, 0, sizeof(*c));
+}
+
+void ut_lstm_free(ut_lstm* l) {
+  ut_free_all(l->W_ih, l->W_hh, l->bias);
+  l->W_ih = l->W_hh = l->bias = NULL;
+}
+
+// x_seq: [T,B,in]. h0/c0: [B,hidden], or NULL for a zero initial state.
+// cache may be NULL for an inference-only pass (no per-step gate state is
+// kept, same as passing NULL to any other layer's _forward in this file).
+// Returns h_seq [T,B,hidden]; writes the final state into *h_n/*c_n if given
+// (both freshly owned tensors the caller must free).
+ut_tensor* ut_lstm_forward_seq(ut_lstm* l, ut_tensor* x_seq, ut_tensor* h0, ut_tensor* c0,
+                               ut_lstm_seq_cache* cache, ut_tensor** h_n, ut_tensor** c_n) {
+  int T = x_seq->shape.shape[0], B = x_seq->shape.shape[1], IN = x_seq->shape.shape[2];
+  int H = l->hidden;
+  ut_sync_cpu(x_seq);
+
+  ut_tensor* h_seq = ut_alloc(3, (int[]){T, B, H}, UT_CPU);
+  if (cache) {
+    cache->steps = malloc((size_t)T * sizeof(ut_lstm_cache));
+    cache->t = T;
+  }
+
+  ut_tensor* h = h0 ? ut_retain(h0) : ut_alloc(2, (int[]){B, H}, UT_CPU);
+  ut_tensor* c = c0 ? ut_retain(c0) : ut_alloc(2, (int[]){B, H}, UT_CPU);
+
+  for (int t = 0; t < T; t++) {
+    ut_tensor* xt = ut_from_data(2, (int[]){B, IN}, x_seq->data + (size_t)t * B * IN, UT_CPU);
+    ut_tensor *c_next, *h_next =
+        ut_lstm_step(l, xt, h, c, &c_next, cache ? &cache->steps[t] : NULL);
+    memcpy(h_seq->data + (size_t)t * B * H, h_next->data, (size_t)B * H * sizeof(float));
+    ut_free_all(xt, h, c);
+    h = h_next;
+    c = c_next;
+  }
+  if (h_n)
+    *h_n = h;
+  else
+    ut_free(h);
+  if (c_n)
+    *c_n = c;
+  else
+    ut_free(c);
+  return h_seq;
+}
+
+// grad_h_seq: [T,B,hidden], the external gradient contribution to h at each
+// timestep (zero everywhere except where h_t is actually consumed — e.g. only
+// the last slice populated for a many-to-one classifier that reads h_n).
+// Accumulates into dW_ih/dW_hh/db; returns dx_seq [T,B,in].
+ut_tensor* ut_lstm_backward_seq(ut_lstm* l, ut_lstm_seq_cache* cache, ut_tensor* grad_h_seq,
+                                ut_tensor* dW_ih, ut_tensor* dW_hh, ut_tensor* db) {
+  int T = cache->t;
+  int B = cache->steps[0].x->shape.shape[0];
+  int IN = l->in, H = l->hidden;
+  ut_sync_cpu(grad_h_seq);
+
+  ut_tensor* dx_seq = ut_alloc(3, (int[]){T, B, IN}, UT_CPU);
+  ut_tensor* dh_next = ut_alloc(2, (int[]){B, H}, UT_CPU);
+  ut_tensor* dc_next = ut_alloc(2, (int[]){B, H}, UT_CPU);
+
+  for (int t = T - 1; t >= 0; t--) {
+    ut_tensor* dh = ut_alloc(2, (int[]){B, H}, UT_CPU);
+    for (int k = 0; k < B * H; k++)
+      dh->data[k] = grad_h_seq->data[(size_t)t * B * H + k] + dh_next->data[k];
+
+    ut_tensor *dx, *dh_prev, *dc_prev;
+    ut_lstm_backward(l, &cache->steps[t], dh, dc_next, dW_ih, dW_hh, db, &dx, &dh_prev, &dc_prev);
+
+    memcpy(dx_seq->data + (size_t)t * B * IN, dx->data, (size_t)B * IN * sizeof(float));
+    ut_free_all(dh, dx, dh_next, dc_next);
+    dh_next = dh_prev;
+    dc_next = dc_prev;
+  }
+  ut_free_all(dh_next, dc_next);
+  return dx_seq;
+}
+
+void ut_lstm_seq_cache_free(ut_lstm_seq_cache* c) {
+  for (int t = 0; t < c->t; t++) ut_lstm_cache_free(&c->steps[t]);
+  free(c->steps);
+  c->steps = NULL;
+  c->t = 0;
+}
+
 // =========================================================
 // Conv2D
 // =========================================================