Project a matrix to the doubly stochastic matrices (Sinkhorn's algorithm)
function B = doubly_stochastic(A)
function B = doubly_stochastic(A, maxiter)
function B = doubly_stochastic(A, [], mode)
function B = doubly_stochastic(A, maxiter, mode)
function B = doubly_stochastic(A, maxiter, mode, checkperiod)
Given an element-wise non-negative matrix A of size nxn, returns a
doubly-stochastic 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^2.
mode (optional):
Setting mode = 1 changes the behavior of the algorithm
such that the input A is an n x p matrix with AA' having
element-wise non-negative entries and the output B is also n x p
such that BB' is a doubly-stochastic matrix. The default value is 0.
checkperiod (optional):
Only check stopping criteria every checkperiod iterations,
to reduce computational burden.0001 function B = doubly_stochastic(A, maxiter, mode, checkperiod) 0002 % Project a matrix to the doubly stochastic matrices (Sinkhorn's algorithm) 0003 % 0004 % function B = doubly_stochastic(A) 0005 % function B = doubly_stochastic(A, maxiter) 0006 % function B = doubly_stochastic(A, [], mode) 0007 % function B = doubly_stochastic(A, maxiter, mode) 0008 % function B = doubly_stochastic(A, maxiter, mode, checkperiod) 0009 % 0010 % Given an element-wise non-negative matrix A of size nxn, returns a 0011 % doubly-stochastic matrix B of size nxn by applying Sinkhorn's algorithm 0012 % to A. 0013 % 0014 % maxiter (optional): 0015 % Strictly positive integer representing the maximum 0016 % number of iterations of Sinkhorn's algorithm. 0017 % The default value of maxiter is n^2. 0018 % mode (optional): 0019 % Setting mode = 1 changes the behavior of the algorithm 0020 % such that the input A is an n x p matrix with AA' having 0021 % element-wise non-negative entries and the output B is also n x p 0022 % such that BB' is a doubly-stochastic matrix. The default value is 0. 0023 % checkperiod (optional): 0024 % Only check stopping criteria every checkperiod iterations, 0025 % to reduce computational burden. 0026 0027 % The file is based on developments in the research paper 0028 % Philip A. Knight, "The Sinkhorn–Knopp Algorithm: Convergence and 0029 % Applications" in SIAM Journal on Matrix Analysis and Applications 30(1), 0030 % 261-275, 2008. 0031 % 0032 % Please cite the Manopt paper as well as the research paper. 0033 0034 % This file is part of Manopt: www.manopt.org. 0035 % Original author: David Young, September 10, 2015. 0036 % Contributors: Ahmed Douik, March 15, 2018. 0037 % Pratik Jawanpuria and Bamdev Mishra, Sep 10, 2019. 0038 % Change log: 0039 % Sep. 10, 2019 (PJ, BM) 0040 % Added the checkperiod parameter. 0041 0042 n = size(A, 1); 0043 tol = eps(n); 0044 0045 if ~exist('maxiter', 'var') || isempty(maxiter) 0046 maxiter = n^2; 0047 end 0048 0049 if ~exist('mode', 'var') || isempty(mode) 0050 mode = 0; 0051 end 0052 0053 if ~exist('checkperiod', 'var') || isempty(checkperiod) 0054 checkperiod = 100; 0055 end 0056 0057 if mode == 1 0058 C = A*A'; 0059 else 0060 C = A; 0061 end 0062 0063 iter = 0; 0064 d_1 = 1./sum(C); 0065 d_2 = 1./(C * d_1.'); 0066 while iter < maxiter 0067 iter = iter + 1; 0068 row = d_2.' * C; 0069 % Check gap condition only at checkperiod intervals. 0070 % It saves computations for large-scale scenarios. 0071 if mod(iter, checkperiod) == 0 0072 gap = max(abs(row .* d_1 - 1)); 0073 if isnan(gap) 0074 break; 0075 end 0076 if gap <= tol 0077 break; 0078 end 0079 end 0080 d_1_prev = d_1; 0081 d_2_prev = d_2; 0082 0083 d_1 = 1./row; 0084 d_2 = 1./(C * d_1.'); 0085 0086 if any(isinf(d_2)) || any(isnan(d_2)) || any(isinf(d_1)) || any(isnan(d_1)) 0087 warning('DoublyStochasticProjection: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 if mode == 1 0096 v = sqrt(d_2(1)/d_1(1))*d_1; 0097 B = diag(v)*A; 0098 else 0099 B = C .* (d_2 * d_1); 0100 end 0101 0102 end