doubly_stochastic_general

 Project a matrix to the doubly stochastic matrices (Sinkhorn's algorithm)

 function [B, d_2, d_1] = doubly_stochastic_general(A, p, q)
 function [B, d_2, d_1] = doubly_stochastic_general(A, p, q, maxiter)
 function [B, d_2, d_1] = doubly_stochastic_general(A, p, q, [], checkperiod)
 function [B, d_2, d_1] = doubly_stochastic_general(A, p, q, maxiter, checkperiod)

 Given an element-wise non-negative matrix A of size nxm, returns a
 matrix B of size nxn by applying Sinkhorn's algorithm
 to A.
 
 maxiter (optional):
    Strictly positive integer representing the maximum 
    number of iterations of Sinkhorn's algorithm. 
    The default value of maxiter is n*m.
 checkperiod (optional):
    Only check stopping criteria every checkperiod iterations,
    to reduce computational burden.

 The file is based on developments in the research paper
 Philip A. Knight, "The Sinkhorn–Knopp Algorithm: Convergence and 
 Applications" in SIAM Journal on Matrix Analysis and Applications 30(1), 
 261-275, 2008.

 Please cite the Manopt paper as well as the above research paper.
 

 See also doubly_stochastic

0001 function [B, d_2, d_1] = doubly_stochastic_general(A, p, q, maxiter, checkperiod)
0002 % Project a matrix to the doubly stochastic matrices (Sinkhorn's algorithm)
0003 %
0004 % function [B, d_2, d_1] = doubly_stochastic_general(A, p, q)
0005 % function [B, d_2, d_1] = doubly_stochastic_general(A, p, q, maxiter)
0006 % function [B, d_2, d_1] = doubly_stochastic_general(A, p, q, [], checkperiod)
0007 % function [B, d_2, d_1] = doubly_stochastic_general(A, p, q, maxiter, checkperiod)
0008 %
0009 % Given an element-wise non-negative matrix A of size nxm, returns a
0010 % matrix B of size nxn by applying Sinkhorn's algorithm
0011 % to A.
0012 %
0013 % maxiter (optional):
0014 %    Strictly positive integer representing the maximum
0015 %    number of iterations of Sinkhorn's algorithm.
0016 %    The default value of maxiter is n*m.
0017 % checkperiod (optional):
0018 %    Only check stopping criteria every checkperiod iterations,
0019 %    to reduce computational burden.
0020 %
0021 % The file is based on developments in the research paper
0022 % Philip A. Knight, "The Sinkhorn–Knopp Algorithm: Convergence and
0023 % Applications" in SIAM Journal on Matrix Analysis and Applications 30(1),
0024 % 261-275, 2008.
0025 %
0026 % Please cite the Manopt paper as well as the above research paper.
0027 %
0028 %
0029 % See also doubly_stochastic
0030 
0031 % This file is part of Manopt: www.manopt.org.
0032 % Original author: Bamdev Mishra, Oct 30, 2020.
0033 % Contributors:
0034 % Change log:
0035 
0036     tol = eps;
0037 
0038     n = size(A, 1);
0039     m = size(A, 2);
0040     
0041     if ~exist('p', 'var') || isempty(p)
0042         p = (1/n)*ones(n,1);
0043     end
0044 
0045     if ~exist('q', 'var') || isempty(q)
0046         q = (1/m)*ones(m,1);
0047     end
0048 
0049     if ~exist('maxiter', 'var') || isempty(maxiter)
0050         maxiter = n*m;
0051     end
0052     
0053     if ~exist('checkperiod', 'var') || isempty(checkperiod)
0054         checkperiod = 100;
0055     end
0056         
0057     C = A;
0058     
0059 
0060     iter = 0;
0061     d_1 = (q')./sum(C); % row vector of size m
0062     d_2 = p./(C * d_1.'); % column vector of size n
0063     gap = inf;
0064     while iter < maxiter
0065         iter = iter + 1;
0066         row = d_2.' * C;
0067 
0068         % Check gap condition only at checkperiod intervals.
0069         % It saves computations for large-scale scenarios.
0070         if mod(iter, checkperiod) == 0 
0071             gap = max(abs(row .* d_1 - q'));
0072             if isnan(gap)
0073                 break;
0074             end
0075             if gap <= tol
0076                 break;
0077             end
0078         end
0079         
0080         d_1_prev = d_1;
0081         d_2_prev = d_2;
0082 
0083         d_1 = (q')./row;
0084         d_2 = p./(C * d_1.');
0085 
0086         if any(isinf(d_2)) || any(isnan(d_2)) || any(isinf(d_1)) || any(isnan(d_1))
0087             warning('DoublyStochasticGeneralProjection:NanInfEncountered', ...
0088                     'Nan or Inf occured at iter %d \n', iter);
0089             d_1 = d_1_prev;
0090             d_2 = d_2_prev;
0091             break;
0092         end
0093     end
0094     
0095     B = C .* (d_2 * d_1);
0096 
0097     fprintf('Iter %d, gap %e \n', iter, gap);
0098 
0099 end