Computationally efficient way to combine degree tables

[tcunis]

The degree table (degmat) of a polynomial describes, for each monomial term, the degree of each variable in that term. Most polynomial operations which manipulate the degree table do so in either of the following ways:

1. Join the degree tables: Operations such as plus or cat join the degree tables of the polynomials involved in the operation, where coefficients, that is, entries in the degree table, that appear in multiple operands are combined (e.g., added or concatenated).
2. Expand and recombine the degree tables: Operations such as times or mtimes expand (permute) the degree tables of the operands to ensure that each coefficient (entry) of any operand is recombined (e.g., multiplied) with all coefficients (entries) of any other operand.

In the current implementation, these use cases involve the following operations:

• For joining, the degree tables are concatenated;
• For expansion and recombination, the degree tables are first permuted via the Kronecker product and added;
• To find duplicate entries in the new degree table, MATLAB's built-in unique function is used; this operation also returns the new degree table in canonical sorting;
• To combine (add) coefficients belonging to duplicate entries, the combined coefficients are multiplied with a matrix of 1 and 0.

While these operations are moderately fast due to their C++ implementation, they are inefficient and computation times become noticeable for larger instances (see, e.g., #63). Let $N_1$ and $N_2$ be the number of degrees in a binary operation; then

• the complexity of the Kronecker product is $\mathcal O(N_1 \cdot N_2)$ (not accounting for the actual size of the coefficients);
• the complexity of unique is $\mathcal O((N_1+N_2) \log(N_1+N_2))$ for 1. and $\mathcal O(N_1 \cdot N_2 \log(N_1 \cdot N_2))$ for 2., assuming an appropriate implementation of the sorting;
• the complexity of the matrix multiplication is $\mathcal O(M \cdot (N_1 + N_2))$ for 1. and $\mathcal O(M \cdot N_1 \cdot N_2)$ for 2., where $M$ is the number of unique degrees, possibly better if sparsity can be exploited.

[tcunis]

Given that the degree tables of the operands are already canonically sorted, joining the tables while restoring the sorting is possible in linear time; expansion and recombination is slightly more expensive.

Here's a pseudo code for combining the degree tables of two polynomials, dga and dgb, along with their coefficients for summation. The number of terms in each polynomial shall be nta and ntb.

ia = ib = 1; % Note: we assume one-based index as in MATLAB
idx = 1; % index to first term of combined degree matrix

degmat = sparse(nta+ntb, nvars); % pre-allocate the degree matrix to the maximal size and number of joined variables
coeffs = casadi.DM(nta+ntb, numel);

while (ia <= nta) || (ib <= ntb)
  % loop over both degree matrices simultaneously
  if (ia > nta) || gt_sort(dga(ia,:), dgb(ib,:)) 
    % only terms of polynomial b left 
    % or the ia-th entry of dga is canonically *after* the ib-th entry of dgb
    degmat(idx,:) = dgb(ib,:);
    coeffs(idx,:) = cfb(ib,:);
    ib++;
  elseif (ib > ntb) || lt_sort(dga(ia,:), dgb(ib,:))
    % only terms of polynomial a left 
    % or the ia-th entry of dga is canonically *before* the ib-th entry of dgb
    degmat(idx,:) = dga(ia,:);
    coeffs(idx,:) = cfa(ia,:);
    ia++;
  else 
    % ia-th entry of dga and ib-th entry of dgb are equal
    degmat(idx,:) = dga(ia,:);
    coeffs(idx,:) = cfa(ia,:) + cfb(ib,:); % replace operation here
    ia++; 
    ib++;
  end

  % advance pointer
  idx++;
end

degmat(idx:end,:) = []; % remove allocated but unused rows of the new degree matrix