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..3f7c0f1 100644 --- a/@myAD/private/matderXmatval.m +++ b/@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 diff --git a/@myAD/private/matvalXmatder.m b/@myAD/private/matvalXmatder.m index d22e859..a52dc29 100644 --- a/@myAD/private/matvalXmatder.m +++ b/@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 diff --git a/compile_mex_files.m b/compile_mex_files.m index 1260377..cef1674 100644 --- a/compile_mex_files.m +++ b/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 ../../; 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