I am not certain that there is any built-in functtion for this sort of problem.
One approach —
M1 = randn(5)
M2 = M1 + randn(5)/10
Mcomp = sqrt((M1-M2).^2/numel(M1));
ms = fitdist(Mcomp(:),'lognormal')
lognparms = @(mu,sigma) [exp(mu + sigma^2/2); exp(2*mu + sigma^2) * (exp(sigma^2)-1); sqrt(exp(2*mu + sigma^2) * (exp(sigma^2)-1))]; % [mean; var; std]
lnparms = lognparms(ms.mu, ms.sigma);
fprintf(1,'\nMean = %.6f\nVar = %.6f\nStDv = %.6f\n',lnparms)
[X,Y] = meshgrid(1:5);
figure
stem3(Mcomp)
hold on
scatter3(X(:), Y(:), Mcomp(:), 50, Mcomp(:), 'filled')
patch([xlim flip(xlim)], [1 1 5 5], ones(1,4)*lnparms(1), 'k', 'FaceAlpha',0.25, 'EdgeColor','none')
hold off
colormap(turbo)
colorbar
axis('padded')
zlabel(["‘Mcomp’ Values For Each" "Comparison Matrix Element"])
The comparison value ‘Mcomp’ is a ‘sort of’ root-mean-square value. (It can be anything, providing it returns a matrix equal to the original matrices.) Since the ‘Mcomp’ calculation I chose can never be negative (and would likely be strictly positive, so non-zero as well), I chose the lognormal distribution to calculate its parameters. If you choose a different calculation, choose an appropriate distribution if you want to calculate its parameters. The gray patch is the mean of the ‘Mcomp’ values, as calculated from the lognormal distribution parameters.
.
