Home > examples > doubly_stochastic_denoising.m

# doubly_stochastic_denoising

## PURPOSE

Find a doubly stochastic matrix closest to a given matrix, in Frobenius norm.

## SYNOPSIS

function doubly_stochastic_denoising()

## DESCRIPTION

``` Find a doubly stochastic matrix closest to a given matrix, in Frobenius norm.

This example demonstrates how to use the geometry factories for the
doubly stochastic multinomial manifold:
multinomialdoublystochasticfactory and
multinomialsymmetricfactory (for the symmetric case.)

The file is based on developments in the research paper
A. Douik and B. Hassibi, "Manifold Optimization Over the Set
of Doubly Stochastic Matrices: A Second-Order Geometry"
ArXiv:1802.02628, 2018.

Please cite the Manopt paper as well as the research paper:
@Techreport{Douik2018Manifold,
Title   = {Manifold Optimization Over the Set of Doubly Stochastic
Matrices: {A} Second-Order Geometry},
Author  = {Douik, A. and Hassibi, B.},
Journal = {Arxiv preprint ArXiv:1802.02628},
Year    = {2018}
}

This can be a starting point for many optimization problems of the form:

minimize f(X) such that X is a doubly stochastic matrix (symmetric or not)

Input:  None. This example file generates random data.

Output: None.

This file is part of Manopt: www.manopt.org.
Original author: Ahmed Douik, March 15, 2018.
Contributors:
Change log:

Xiaowen Jiang Aug. 31, 2021

## CROSS-REFERENCE INFORMATION

This function calls:
This function is called by:

## SOURCE CODE

```0001 function doubly_stochastic_denoising()
0002 % Find a doubly stochastic matrix closest to a given matrix, in Frobenius norm.
0003 %
0004 % This example demonstrates how to use the geometry factories for the
0005 % doubly stochastic multinomial manifold:
0006 %  multinomialdoublystochasticfactory and
0007 %  multinomialsymmetricfactory (for the symmetric case.)
0008 %
0009 % The file is based on developments in the research paper
0010 % A. Douik and B. Hassibi, "Manifold Optimization Over the Set
0011 % of Doubly Stochastic Matrices: A Second-Order Geometry"
0012 % ArXiv:1802.02628, 2018.
0013 %
0014 % Link to the paper: https://arxiv.org/abs/1802.02628.
0015 %
0016 % Please cite the Manopt paper as well as the research paper:
0017 % @Techreport{Douik2018Manifold,
0018 %   Title   = {Manifold Optimization Over the Set of Doubly Stochastic
0019 %              Matrices: {A} Second-Order Geometry},
0020 %   Author  = {Douik, A. and Hassibi, B.},
0021 %   Journal = {Arxiv preprint ArXiv:1802.02628},
0022 %   Year    = {2018}
0023 % }
0024 %
0025 % This can be a starting point for many optimization problems of the form:
0026 %
0027 % minimize f(X) such that X is a doubly stochastic matrix (symmetric or not)
0028 %
0029 % Input:  None. This example file generates random data.
0030 %
0031 % Output: None.
0032 %
0033 % This file is part of Manopt: www.manopt.org.
0034 % Original author: Ahmed Douik, March 15, 2018.
0035 % Contributors:
0036 % Change log:
0037 %
0038 %    Xiaowen Jiang Aug. 31, 2021
0040
0041     %% Generate input data
0042     n = 100;
0043     sigma = 1/n^2;
0044     % Generate a doubly stochastic matrix using the Sinkhorn algorithm
0045     B = doubly_stochastic(abs(randn(n, n)));
0046     % Add noise to the matrix
0047     A = max(B + sigma*randn(n, n), 0.01);
0048
0049     %% Setup an optimization problem for denoising
0050
0051     % Manifold initialization: pick the symmetric or non-symmetric case
0052     symmetric_case = true;
0053     if ~symmetric_case
0054         manifold = multinomialdoublystochasticfactory(n);
0055     else
0056         % Input must also be symmetric (otherwhise, gradient must be adapted.)
0057         A = (A+A')/2;
0058         manifold = multinomialsymmetricfactory(n);
0059     end
0060     problem.M = manifold;
0061
0062     % Specify cost function and derivatives
0063     problem.cost = @(X) 0.5*norm(A-X, 'fro')^2;
0065     problem.ehess = @(X, U) U;   % Euclidean Hessian
0066
0067     % An alternative way to compute the egrad and the ehess is to use
0068     % automatic differentiation provided in the deep learning toolbox (slower).
0069     % Notice that the function norm is not supported for AD so far.
0070     % Replace the expression norm(..., 'fro')^2 with cnormsqfro described
0071     % in the file manoptADhelp.m.
0072     % problem.cost = @(X) 0.5*cnormsqfro(A-X);
0075
0076     % Check consistency of the gradient and the Hessian. Useful if you
0077     % adapt this example for a new cost function and you would like to make
0078     % sure there is no mistake. These checks are optional of course.
0079     warning('off', 'manopt:multinomialdoublystochasticfactory:exp');
0080     warning('off', 'manopt:multinomialsymmetricfactory:exp');
0081     figure();
0083     figure();
0084     checkhessian(problem);  % Check the hessian. This test usually fails
0085                             % due to the use of a first order retraction,
0086                             % unless the test is performed at a critical point.
0087
0088     %% Solve
0089     % A first order algorithm
0090     % Minimize the cost function using the Conjugate Gradients algorithm.
0091     [X1, ~, info, ~] = conjugategradient(problem);
0093
0094     % Computing the distance between the original, noiseless matrix and the
0095     % found solution
0096     fprintf('||X-B||_F^2 = %g\n', norm(X1 - B, 'fro')^2);
0097
0098     % A second order algorithm
0099     % Minimize the cost function using the trust-regions method.
0100     [X2, ~, info, ~] = trustregions(problem);
0102
0103     % Computing the distance between the original, noiseless matrix and the
0104     % found solution
0105     fprintf('||X-B||_F^2 = %g\n', norm(X2 - B, 'fro')^2);
0106
0107
0108     %% Display
0109     figure();
0110     semilogy(S1, 'Color', 'red', 'Marker', '*', 'LineWidth', 2, ...
0112     hold on;
0113     semilogy(S2, 'Color', 'blue', 'Marker', '+', 'LineWidth', 2, ...
0114                  'DisplayName', 'Trust Regions Method');
0115
0116     set(gca, 'YGrid', 'on', 'XGrid', 'on');
0117     xlabel('Number of Iteration', 'FontSize', 16);
0118     ylabel('Norm of Gradient', 'FontSize', 16);
0119     legend1 = legend('show');
0120     set(legend1, 'FontSize', 16);
0121
0122     % This Hessian test is performed at a computed solution, which is
0123     % expected to be a critical point. Thus we expect the slope test to
0124     % succeed. It could fail if the solution X1 has entries very close to
0125     % zero, so that numerical issues might come up. This is because de
0126     % Fisher metric on the multinomial manifold involves division by the
0127     % entries of X.
0128     figure();
0129     checkhessian(problem, X1);
0130
0131 end```

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