From d0c181f7b37d3a7fb1825d3658d04a5b808a6d5e Mon Sep 17 00:00:00 2001 From: gregkaplan Date: Wed, 29 Jul 2026 17:20:16 +1000 Subject: [PATCH 1/3] Speed up sparse kernels: kron-based matrix products, factorization reuse in mldivide, add blkdiag - matvalXmatder / matderXmatval: replace the explicit loop over the inner dimension with a single kron-based sparse product (identities vec(A*dB) = kron(I,A)*vec(dB) and vec(dA*B) = kron(B',I)*vec(dA)). 30-150x faster at kernel level; package-level mtimes for (1000x1000 sparse)*(1000x20) with 100 derivative directions: AD*AD 0.196s -> 0.040s, double*AD 0.125s -> 0.0008s. Results identical. - mldivide: factor the matrix once per call via decomposition() and reuse it for the value solve and all derivative solves. Previously each backslash refactorized: twice for one RHS, once per column for matrix RHS. 13% (1 RHS) to 27% (10 RHS) faster on a 3600x3600 5-point stencil; gains grow with factorization cost. - blkdiag: new method assembling values and stacked derivatives with one sparse() call each. Building block-diagonal AD matrices by repeated concatenation is O(nblocks^2) with a full derivative permutation per step; blkdiag of 30 blocks (40x40, 100 directions) is 68x faster with bit-identical output. - compile_mex_files.m: comment out remaining mex compilation with an explanatory note. On recent MATLAB (tested R2025b) the .m fallbacks are as fast or faster (matdrivXvecval mex measured 4-37x slower than its .m), and a compiled mex silently shadows the .m implementation. All changed paths validated against central finite differences (~1e-10) and on a continuous-time HANK transition-Jacobian application (198x198 Jacobian through backward HJB / forward KFE with sparse solves): Jacobians identical to master to the last bit. Co-Authored-By: Claude Fable 5 --- @myAD/blkdiag.m | 61 +++++++++++++++++++++++++++++++++++ @myAD/mldivide.m | 31 +++++++++++------- @myAD/private/matderXmatval.m | 50 +++++++++------------------- @myAD/private/matvalXmatder.m | 51 +++++++++-------------------- compile_mex_files.m | 20 +++++++++--- 5 files changed, 127 insertions(+), 86 deletions(-) create mode 100644 @myAD/blkdiag.m diff --git a/@myAD/blkdiag.m b/@myAD/blkdiag.m new file mode 100644 index 0000000..9c21a4c --- /dev/null +++ b/@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 diff --git a/@myAD/mldivide.m b/@myAD/mldivide.m index 73653af..288f7cd 100644 --- a/@myAD/mldivide.m +++ b/@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; diff --git a/@myAD/private/matderXmatval.m b/@myAD/private/matderXmatval.m index 88355c9..a5f42de 100644 --- a/@myAD/private/matderXmatval.m +++ b/@myAD/private/matderXmatval.m @@ -1,40 +1,20 @@ 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) + % + % Uses the identity vec(dA*B) = kron(B', I_nrow)*vec(dA), so the whole + % operation is a single sparse matrix product. This replaces an explicit + % loop over the inner dimension which was O(ninter*(nnz(A)+nnz(B))) and + % dominated runtime for large matrices. + % % by SeHyoun Ahn, July 2018 + % updated to kron formulation, 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); + nrow = size(A, 1)/size(B, 1); + output = kron(sparse(B'), speye(nrow))*A; end diff --git a/@myAD/private/matvalXmatder.m b/@myAD/private/matvalXmatder.m index d22e859..92c2117 100644 --- a/@myAD/private/matvalXmatder.m +++ b/@myAD/private/matvalXmatder.m @@ -1,40 +1,21 @@ 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) + % + % Uses the identity vec(A*dB) = kron(I_ncol, A)*vec(dB), so the whole + % operation is a single sparse matrix product. This replaces an explicit + % loop over the inner dimension which was O(ninter*(nnz(A)+nnz(B))) and + % dominated runtime for large matrices. + % % by SeHyoun Ahn, July 2018 + % updated to kron formulation, July 2026 - [Arow, Acol, Aval] = find(A); - [nrow, ninter] = size(A); - Arow = Arow(:)'; - Acol = Acol(:)'; - Aval = Aval(:)'; - - [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); + ninter = size(A, 2); + ncol = size(B, 1)/ninter; + output = kron(speye(ncol), sparse(A))*B; end diff --git a/compile_mex_files.m b/compile_mex_files.m index 1260377..59813b6 100644 --- a/compile_mex_files.m +++ b/compile_mex_files.m @@ -1,6 +1,18 @@ -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; +% NOTE (July 2026): On recent MATLAB releases (tested R2025b), the pure-MATLAB +% implementations in @myAD/private are as fast as or faster than these mex +% kernels, which predate MATLAB's multithreaded sparse operations and implicit +% expansion: +% - valXder.m is a one-line implicit-expansion product (parity with mex) +% - matdrivXvecval.m was benchmarked 4x-37x FASTER than the compiled mex +% - matvalXmatder.m / matderXmatval.m use a single kron-based sparse product +% Since a compiled mex file shadows the .m file of the same name, compiling +% these can make the package slower. Compilation is left here only for use on +% old MATLAB releases; delete stale compiled binaries (*.mexa64, *.mexmaci64, +% *.mexmaca64, *.mexw64) from @myAD/private to use the .m implementations. +% +% 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; % mex -v -largeArrayDims COMPFLAGS='-O2 -ftree-vectorize -ftree-vectorize-verbose=7 -fopt-info-missed -Wall' dertransp.c; -cd ../../; +% cd ../../; From 7150b1ec1811bb163838024976beee1f62876958 Mon Sep 17 00:00:00 2001 From: gregkaplan Date: Thu, 30 Jul 2026 10:58:41 +1000 Subject: [PATCH 2/3] Make matrix kernels memory-safe and complex-correct; add regression tests Follow-up to the previous commit, addressing two problems with the kron-based formulation. Memory: kron(I,A) and kron(B',I) materialize ncol*nnz(A) and nrow*nnz(B) entries respectively, regardless of how sparse the derivative payload is. For a dense 2000x2000 value matrix with 100 columns that is 4e8 stored entries (~5 GB), even if the payload has a handful of nonzeros. Both kernels now avoid forming any kron intermediate: - matvalXmatder: the contraction is over the leading index of the column-wise stacking, so reshaping the payload to (ninter x ncol*nderiv) lays every direction side by side and one sparse product handles all of them. Strictly better than kron: same result, O(nnz) peak memory, and faster everywhere measured (dense-A/sparse-payload case 0.0067s -> 0.0001s; HANK-like case 0.0007s -> 0.0003s). - matderXmatval: the contraction is over the trailing index, so the payload is reindexed once (O(nnz)) to move the derivative direction into the row index, one sparse product contracts all directions, and the result is reindexed back. Faster than kron in the wide/dense regimes (5-6x) and within ~30% of it in the moderate-sparse regime, with bounded memory. Complex correctness: neither kernel now takes a transpose at all, so both are correct for complex input. The kron version of matderXmatval used B' (conjugate transpose) where vec(dA*B) = kron(B.',I)*vec(dA) requires the nonconjugate transpose. Note master had an analogous pre-existing bug in both kernels: `Aval = Aval(:)'` conjugates the stored derivative values. Measured max derivative errors for complex input, versus a per-direction reference: master 5.1 and 7.3, kron-with-B' 4.2, both new kernels ~5e-16. Adds test_matrix_kernels.m: checks mtimes (AD*AD, AD*double, double*AD) and mldivide against the per-direction product rule over five shapes, for real and complex inputs, plus zero-payload, single-direction, and blkdiag cases. The complex cases fail on master, which is how the conjugation bug was found. Co-Authored-By: Claude Fable 5 --- @myAD/private/matderXmatval.m | 42 ++++++++-- @myAD/private/matvalXmatder.m | 24 ++++-- test_matrix_kernels.m | 152 ++++++++++++++++++++++++++++++++++ 3 files changed, 204 insertions(+), 14 deletions(-) create mode 100644 test_matrix_kernels.m diff --git a/@myAD/private/matderXmatval.m b/@myAD/private/matderXmatval.m index a5f42de..3f7c0f1 100644 --- a/@myAD/private/matderXmatval.m +++ b/@myAD/private/matderXmatval.m @@ -7,14 +7,42 @@ % % Output: derivative of A*B stacked column-wise: (nrow*ncol x nderiv) % - % Uses the identity vec(dA*B) = kron(B', I_nrow)*vec(dA), so the whole - % operation is a single sparse matrix product. This replaces an explicit - % loop over the inner dimension which was O(ninter*(nnz(A)+nnz(B))) and - % dominated runtime for large matrices. + % 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 - % updated to kron formulation, July 2026 + % reformulated July 2026 + + [ninter, ncol] = size(B); + 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); - nrow = size(A, 1)/size(B, 1); - output = kron(sparse(B'), speye(nrow))*A; + % 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 diff --git a/@myAD/private/matvalXmatder.m b/@myAD/private/matvalXmatder.m index 92c2117..a52dc29 100644 --- a/@myAD/private/matvalXmatder.m +++ b/@myAD/private/matvalXmatder.m @@ -7,15 +7,25 @@ % % Output: derivative of A*B stacked column-wise: (nrow*ncol x nderiv) % - % Uses the identity vec(A*dB) = kron(I_ncol, A)*vec(dB), so the whole - % operation is a single sparse matrix product. This replaces an explicit - % loop over the inner dimension which was O(ninter*(nnz(A)+nnz(B))) and - % dominated runtime for large matrices. + % 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 - % updated to kron formulation, July 2026 + % reformulated July 2026 - ninter = size(A, 2); + [nrow, ninter] = size(A); + nderiv = size(B, 2); ncol = size(B, 1)/ninter; - output = kron(speye(ncol), sparse(A))*B; + + output = reshape(sparse(A)*reshape(B, ninter, ncol*nderiv), nrow*ncol, nderiv); end diff --git a/test_matrix_kernels.m b/test_matrix_kernels.m new file mode 100644 index 0000000..7196d70 --- /dev/null +++ b/test_matrix_kernels.m @@ -0,0 +1,152 @@ +% Regression tests for the matrix-product derivative kernels reached through +% mtimes and mldivide (private/matvalXmatder, private/matderXmatval) and for +% blkdiag. +% +% For R = X*Y the derivative in direction k must satisfy +% unstack(dR(:,k)) = unstack(dX(:,k))*Yval + Xval*unstack(dY(:,k)), +% where unstack reshapes a column-wise stacked column back into a matrix. Each +% test builds myAD objects with explicit derivative matrices and compares against +% that definition, evaluated one direction at a time. +% +% Complex cases are included because the vectorization identities behind the +% kernels require nonconjugate transposes: a conjugate transpose (or a conjugating +% ctranspose applied to the stored values) passes every real-valued test and +% silently returns wrong derivatives for complex input. +% +% July 2026 + +clear; +warning('off', 'AutoDiff:autoext'); +rng(20260730); + +tol = 1e-11; +n_fail = 0; +fprintf('%-42s %10s %s\n', 'case', 'max err', 'status'); + +shapes = {[4 3 2], [7 5 6], [1 4 3], [5 1 2], [6 6 1]}; + +for scenario = {'real', 'complex'} + is_complex = strcmp(scenario{1}, 'complex'); + + for s = 1:numel(shapes) + n = shapes{s}(1); m = shapes{s}(2); p = shapes{s}(3); + K = 5; + + Xval = randn(n,m); Yval = randn(m,p); + dX = sprandn(n*m, K, 0.7); dY = sprandn(m*p, K, 0.7); + if is_complex + Xval = Xval + 1i*randn(n,m); Yval = Yval + 1i*randn(m,p); + dX = dX + 1i*sprandn(n*m, K, 0.7); dY = dY + 1i*sprandn(m*p, K, 0.7); + end + + X = myAD(Xval, dX); + Y = myAD(Yval, dY); + + ref_both = zeros(n*p, K); + ref_xad = zeros(n*p, K); + ref_yad = zeros(n*p, K); + for k = 1:K + dXk = reshape(full(dX(:,k)), n, m); + dYk = reshape(full(dY(:,k)), m, p); + ref_both(:,k) = reshape(dXk*Yval + Xval*dYk, [], 1); + ref_xad(:,k) = reshape(dXk*Yval, [], 1); + ref_yad(:,k) = reshape(Xval*dYk, [], 1); + end + + got_both = getderivs(X*Y); + got_xad = getderivs(X*Yval); + got_yad = getderivs(Xval*Y); + + [n_fail] = report(n_fail, sprintf('mtimes AD*AD %s %dx%dx%d', scenario{1}, n, m, p), got_both, ref_both, tol); + [n_fail] = report(n_fail, sprintf('mtimes AD*dbl %s %dx%dx%d', scenario{1}, n, m, p), got_xad, ref_xad, tol); + [n_fail] = report(n_fail, sprintf('mtimes dbl*AD %s %dx%dx%d', scenario{1}, n, m, p), got_yad, ref_yad, tol); + end + + % ---- mldivide with matrix right-hand side, both AD ---- + nq = 5; pq = 3; K = 4; + Aval = randn(nq,nq) + 4*eye(nq); Bval = randn(nq,pq); + dA = sprandn(nq*nq, K, 0.7); dB = sprandn(nq*pq, K, 0.7); + if is_complex + Aval = Aval + 1i*randn(nq,nq); Bval = Bval + 1i*randn(nq,pq); + dA = dA + 1i*sprandn(nq*nq, K, 0.7); dB = dB + 1i*sprandn(nq*pq, K, 0.7); + end + A = myAD(Aval, dA); B = myAD(Bval, dB); + Zval = Aval\Bval; + ref = zeros(nq*pq, K); + for k = 1:K + dAk = reshape(full(dA(:,k)), nq, nq); + dBk = reshape(full(dB(:,k)), nq, pq); + ref(:,k) = reshape(Aval\(dBk - dAk*Zval), [], 1); + end + got = getderivs(A\B); + [n_fail] = report(n_fail, sprintf('mldivide AD\\AD %s', scenario{1}), got, ref, 1e-9); +end + +% ---- edge case: all-zero derivative payload, shape preserved ---- +n = 4; m = 3; p = 2; +X = myAD(randn(n,m), sparse(n*m, 3)); +Y = myAD(randn(m,p), sparse(m*p, 3)); +got = getderivs(X*Y); +err = full(max(abs(got(:)))); +ok = (err == 0) && isequal(size(got), [n*p, 3]); +n_fail = n_fail + ~ok; +fprintf('%-42s %10.2e %s\n', 'zero payload, shape preserved', err, status_str(ok)); + +% ---- edge case: a single derivative direction ---- +Xval = randn(n,m); Yval = randn(m,p); +dX = sprandn(n*m, 1, 0.8); +X = myAD(Xval, dX); +ref = reshape(reshape(full(dX), n, m)*Yval, [], 1); +[n_fail] = report(n_fail, 'single derivative direction', getderivs(X*Yval), ref, tol); + +% ---- blkdiag against explicit assembly ---- +K = 6; +b1 = myAD(randn(3,2), sprandn(6, K, 0.6)); +d2 = randn(2,2); +b3 = myAD(randn(4,3), sprandn(12, K, 0.6)); +Ablk = blkdiag(b1, d2, b3); + +Nr = 9; Nc = 7; +val_ref = zeros(Nr, Nc); +val_ref(1:3, 1:2) = getvalues(b1); +val_ref(4:5, 3:4) = d2; +val_ref(6:9, 5:7) = getvalues(b3); +err_v = full(max(abs(getvalues(Ablk) - val_ref), [], 'all')); + +der_ref = zeros(Nr*Nc, K); +spec = {getderivs(b1), 0, 0, 3, 2; getderivs(b3), 5, 4, 4, 3}; +for bi = 1:size(spec,1) + der = spec{bi,1}; r0 = spec{bi,2}; c0 = spec{bi,3}; rb = spec{bi,4}; cb = spec{bi,5}; + for k = 1:K + full_k = zeros(Nr, Nc); + full_k(r0+(1:rb), c0+(1:cb)) = reshape(full(der(:,k)), rb, cb); + der_ref(:,k) = der_ref(:,k) + reshape(full_k, [], 1); + end +end +err_d = full(max(abs(getderivs(Ablk) - der_ref), [], 'all')); +ok = (err_v == 0) && (err_d < tol); +n_fail = n_fail + ~ok; +fprintf('%-42s %10.2e %s\n', 'blkdiag vs explicit assembly', max(err_v, err_d), status_str(ok)); + +fprintf('\n'); +if n_fail == 0 + fprintf('All matrix-kernel tests passed.\n'); +else + error('%d matrix-kernel test(s) FAILED.', n_fail); +end + + +function n_fail = report(n_fail, name, got, ref, tol) + err = full(max(abs(got(:) - ref(:)))); + ok = err < tol; + n_fail = n_fail + ~ok; + fprintf('%-42s %10.2e %s\n', name, err, status_str(ok)); +end + +function s = status_str(ok) + if ok + s = 'ok'; + else + s = 'FAIL'; + end +end From e161ed2d3449af351294f90682f0162a38884d59 Mon Sep 17 00:00:00 2001 From: gregkaplan Date: Thu, 30 Jul 2026 13:32:21 +1000 Subject: [PATCH 3/3] Correct the compile_mex_files note: matvalXmatder.c uses an older calling convention The previous note said the mex kernels were merely slower, and described the new matvalXmatder/matderXmatval as kron-based (no longer true after the reformulation). The more important point was missing: matvalXmatder.c expects the derivative of the TRANSPOSED right-hand matrix, as its own header documents and as the old mtimes.m supplied via dertransp. The 2018 .m takes the column-wise stacked derivative instead, and today's mtimes.m calls it that way. Because a compiled mex silently shadows the same-named .m, compiling matvalXmatder.c makes AD matrix-matrix products return wrong derivatives with no warning. Verified: with the column-wise convention the mex is off by O(1); called as mex(A, dertransp(dB, n)) it matches the .m exactly. So the kernel is not buggy in itself -- it is correct through the old API -- but it must not be compiled against the current mtimes. Co-Authored-By: Claude Fable 5 --- compile_mex_files.m | 40 ++++++++++++++++++++++++++++------------ 1 file changed, 28 insertions(+), 12 deletions(-) diff --git a/compile_mex_files.m b/compile_mex_files.m index 59813b6..cef1674 100644 --- a/compile_mex_files.m +++ b/compile_mex_files.m @@ -1,18 +1,34 @@ -% NOTE (July 2026): On recent MATLAB releases (tested R2025b), the pure-MATLAB -% implementations in @myAD/private are as fast as or faster than these mex -% kernels, which predate MATLAB's multithreaded sparse operations and implicit -% expansion: -% - valXder.m is a one-line implicit-expansion product (parity with mex) -% - matdrivXvecval.m was benchmarked 4x-37x FASTER than the compiled mex -% - matvalXmatder.m / matderXmatval.m use a single kron-based sparse product -% Since a compiled mex file shadows the .m file of the same name, compiling -% these can make the package slower. Compilation is left here only for use on -% old MATLAB releases; delete stale compiled binaries (*.mexa64, *.mexmaci64, -% *.mexmaca64, *.mexw64) from @myAD/private to use the .m implementations. +% 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; +% 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 ../../;