Home > manopt > manifolds > multinomial > doubly_stochastic_general.m

# doubly_stochastic_general

## PURPOSE

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

## SYNOPSIS

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

## DESCRIPTION

``` 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```

## CROSS-REFERENCE INFORMATION

This function calls:
This function is called by:

## SOURCE CODE

```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```

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005