Home > manopt > manifolds > multinomial > doubly_stochastic.m

# doubly_stochastic

## PURPOSE

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

## SYNOPSIS

function B = doubly_stochastic(A, maxiter, mode, checkperiod)

## DESCRIPTION

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

## CROSS-REFERENCE INFORMATION

This function calls:
This function is called by:

## SOURCE CODE

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

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