Asked by Giuseppe Pintori
on 13 Sep 2019 at 10:30

Hi guys, I need your help.

I want to create a matrix(4,4) in which the main diagonal have values between 0.3 and 1 and the other cells assume values such as to have a horizontal sum equal to 1.

By now I'm using the following code but the only result is to have a main diagonal composed by the same numbers:

x = eye(4)

x(1,1) = 1+(0.3-1)*rand(1,1)

x(2,2) = x(1,1)

x(3,3) = x(1,1)

x(4,4) = x(1,1)

Any suggestion?

PS : I've tried even with diag

Answer by Stephen Cobeldick
on 13 Sep 2019 at 11:08

Edited by Stephen Cobeldick
on 13 Sep 2019 at 12:49

Accepted Answer

>> N = 4; % matrix size

>> M = nan(N,N); % preallocate

>> V = 0.3+(1-0.3)*rand(1,1) % diagonal value

V =

0.47505

>> M(~eye(N)) = randfixedsum(N-1,N,1-V,0,1); % other values

>> M = M.'; % transpose

>> M(1:N+1:end) = V % assign diagonal value

M =

0.47505 0.40657 0.0087969 0.10958

0.12917 0.47505 0.21287 0.1829

0.35794 0.15519 0.47505 0.011825

0.41335 0.032696 0.078907 0.47505

Checking the sum of each row:

>> sum(M,2)

ans =

1

1

1

1

and diagonal:

>> diag(M)

ans =

0.47505

0.47505

0.47505

0.47505

Bruno Luong
on 13 Sep 2019 at 14:32

@John, I have some formula on the distribution for you (you might know it already).

for X, m vectors of R^n as returned by

X = randfixedsum(n,m,1,0,1)

X(:,i) = { x uniform in [0,1]^n such that sum(x) = 1 }, i=1...,m.

the projection of any component #k, x(k,:) is on [0,1] has the pdf that is a polynomial of order (n-2):

P(x) := (n-1)(1-x).^(n-2) = sum_ {k=0,1..., n-2} p_{n,k} x^k;

with p_{n,k} = (n-1)*(-1)^k* nchoosek(n-2,k);

Example:

n=4;

m = 1e6;

X = randfixedsum(n,m,1,0,1);

P=arrayfun(@(k) (n-1)*(-1)^k*nchoosek(n-2,k),n-2:-1:0);

figure()

ax1 = subplot(2,1,1);

ezplot(ax1,@(x) (n-1)*(1-x).^(n-2),[0 1]);

legend(ax1,"PDF n="+n)

ax2 = subplot(2,1,2);

histogram(X(:),100); % by symmetry all components have the same PDF

linkaxes([ax1 ax2],'x');

xlim(ax1,[0 1]);

ylim(ax1,[0 n-1]);

The brutefoce rejection method from Stephen's (before-edited) code implies the diagonal elements follow the PDF of polynomial form of 2nd order 3*(1-x)^2, truncated bellow 0.3.

In your case as you generate the diagonal as uniform, therefore it's certainly different than the distribution from Stephen's method (polynomial).

The PDF formula can be showed recursively on n, without much difficulty technically.

John D'Errico
on 13 Sep 2019 at 17:43

Yes. I had come to that conclusion after some thought too. But there still seems to be the question of what is the true goal, since only Giuseppe can know that. I think Stephen's solution comes closer than mine.

On the other hand, HAD I generated the diagonal elements using a better distribution than uniform, then my solution would be an avenue to not needing to use a rejection while loop at all.

Bruno Luong
on 13 Sep 2019 at 19:39

Oh I see your point. I think you are right by forcing the right PDF to generate D, it's equivalent to rejection method. This I'm pretty sure it's right because of the special property of simplex and barycentric coordinates.

So the right way (I simplify the procedure to a single row without loosing the generality) is

N = 4;

dmin = 0.3;

dmax = 1;

d = dmax-(dmax-dmin)*rand().^(1/(N-1)); % rather than dmin+(dmax-dmin)*rand;

v = (1-d) * randfixedsum(N-1,1,1,0,1);

then insert d to v....

Sign in to comment.

Answer by John D'Errico
on 13 Sep 2019 at 11:19

Edited by John D'Errico
on 13 Sep 2019 at 11:21

Easy enough, it seems. First, determine the diagonal elements.

x = diag(rand(1,4)*.7 + .3);

Next, you need to choose the other row elements randomly so the sum will be 1. But that sum will now depend on the diagonal element you just chose. Stilll simple, as long as you use randfixedsum, by Roger Stafford, found on the file exchange.

for i = 1:4

x(i,setdiff(1:4,i)) = randfixedsum(3,1,1 - x(i,i),0,1)';

end

Did it work? Of course.

x

x =

0.83586 0.075979 0.057706 0.030454

0.012356 0.85664 0.11425 0.016757

0.13748 0.21163 0.43081 0.22009

0.15838 0.037488 0.16129 0.64284

>> sum(x,2)

ans =

1

1

1

1

Find randfixedsum here:

Sign in to comment.

Answer by Matt J
on 13 Sep 2019 at 10:42

x=eye(4);

x(1:5:end)=0.7*rand(4,1)+0.3

Giuseppe Pintori
on 13 Sep 2019 at 12:18

sorry, I forgot to write that the diagonal must be composed of identical values

Sign in to comment.

Answer by Bruno Luong
on 13 Sep 2019 at 17:36

Edited by Bruno Luong
on 13 Sep 2019 at 20:04

Here is a method that has two advantages:

- without the need of Roger's FEX randfixedsum
- Produce matrix with rigourous uniform conditional probability

N = 4; % matrix size

% diagonal lo/up bounds

dmin = 0.3;

dmax = 1;

% random (common) diagonal value

d = dmax-(dmax-dmin)*rand().^(1/(N-1)); % Edit see comment above, equiv to rejection method

% d = dmin+(dmax-dmin)*rand;

% Generate N random vectors of length N-1 required sum == (1-d)

V = -log(rand(N-1,N)); % Marsaglia's [1961] method

V = V .* ((1-d)./sum(V,1));

% Arrange in the final matrix

A = zeros(N);

isdiag = sparse(1:N,1:N,true);

A(isdiag) = d;

A(~isdiag) = V(:);

A = A.';

% Check result

disp(A)

sum(A,2)

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 2 Comments

## Stephen Cobeldick (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/480247-main-diagonal-operations-problem#comment_745526

## Giuseppe Pintori (view profile)

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/480247-main-diagonal-operations-problem#comment_745571

Sign in to comment.