Instantly share code, notes, and snippets.

Embed
What would you like to do?
Matlab implementation of the "deleted pivot" operation described by Greg Kuperberg in "Kasteleyn cokernels"
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment