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.
Link to the paper: https://arxiv.org/abs/1802.02628.
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
Added AD to compute the egrad and the ehess0001 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 0039 % Added AD to compute the egrad and the ehess 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; 0064 problem.egrad = @(X) X - A; % Euclidean gradient 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); 0073 % Call manoptAD to prepare AD for the problem structure 0074 % problem = manoptAD(problem); 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(); 0082 checkgradient(problem); % Check the gradient 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); 0092 S1 = [info.gradnorm] ; % Collecting the Gradient information 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); 0101 S2 = [info.gradnorm] ; % Collecting the Gradient information 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, ... 0111 'DisplayName', 'Conjugate Gradient Algorithm'); 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