Matlab implementation of the "deleted pivot" operation described by Greg Kuperberg in "Kasteleyn cokernels"

matrix_deleted_pivot.m

1;  % script file

function res = deleted_pivot(M, i, j)
    % Calculates the deleted pivot of matrix M at row i, column j
    % as described by Greg Kuperberg in "Kasteleyn cokernels"

    % X = the jth column of M with element i removed
    X = [M(:,j)(1:i-1); M(:,j)(i+1:end)];
    % Yt = the ith row of M with element j removed
    Yt = [M(i,:)(1:j-1), M(i,:)(j+1:end)];
    % Mxy is what we're calling M', which is the matrix M with row i
    % and column j removed.
    Mx = [M(:, 1:j-1), M(:, j+1:end)];
    Mxy = [Mx(1:i-1, :); Mx(i+1:end, :)];
    % The result is given as M' - X*Yt, where * is the outer product
    res = Mxy - X*Yt;
endfunction

function res = deleted_pivots(M, varargin)
    % Calculates the deleted pivot of matrix M at rows i_k, and columns
    % j_k given by arguments [i_k, j_k].
    % Example:
    %   M = eye(5)
    %   Mdp = deleted_pivots(M, [2, 4], [1, 3], [5, 5])

    % Stop recursion if we're out of varargin
    if (numel(varargin) == 0)
        res = M;
        return;
    endif
    % Validate and then do the actual computation
    ij = varargin{1};
    assert(numel(ij) == 2, "Arguments to deleted_pivots must be two-element arrays.")
    M_new = deleted_pivot(M, ij(1), ij(2));
    % Compute varargin for the next recursion
    next_varargin = varargin(2:end);
    reduce_indices_curried = @(i) @(j) @(indices) reduce_indices(i, j, indices);
    reduce_indices_ij = reduce_indices_curried(ij(1))(ij(2));
    next_varargin = cellfun(reduce_indices_ij, next_varargin, "UniformOutput", false);
    % Recurse!
    res = deleted_pivots(M_new, next_varargin{:});
endfunction

function res = reduce_indices(i, j, indices)
    % Reduce indices [i_k, j_k] by one if they are greater than i or j,
    % respectively. This is useful for updating indices when a row or
    % column of a matrix has been removed.
    res = indices;
    if (indices(1) > i)
        res(1) -= 1;
    endif
    if (indices(2) > j)
        res(2) -= 1;
    endif
endfunction