How to generate a semipositive indecomposible matrix of nxn

[anthonyQueen]

I am still using MC 11.

I was able to make a random square matrix called A of n dimension, with 1<n<1000 (infinite would have been better but my pc horse power is limited). A needed to be positive, with Max Eigenvalue  lower than 1.

I used the following combination of formulas, which probably is not the optimum, but it works.

1): Z(z,z)=:rnd(1/n)    # so to have numbers that fulfill the condition MaxEval<1

2) A=: matrix(n,n,Z)

I would like to refine it, and I ask kindly advice.

The constrains that I would like to add  in generating A are the following.

a) Still a nxn square Matrix with Max Eigenvalue lower than 1.

b) But A should be semi-positive, with random coefficients, including zeros and not just positive as above.

c) However A should be indecomposable (called also irreducible), i.e. A^x  (with x a positive number) should turn a positive matrix.

Thanks for any suggestion.

[LucMeekes]

Use Z (z,z)=:floor (rnd (n))/n^2.
That will cause a 0 every now and then.

Success!
Luc

[anthonyQueen]

Many Thanks Luc!

Seems simple and working. Yet I noticed that it generates only n numbers, zeros included. So for a matrix 3x3 there are only 3 numbers generated, which are put in the 9 cells.  And this may be a problem, since each cell does not have a random semipositive number (I should probably mentioned that in my constraints). I fear that  matrix A could become singular, which is not what I want. 

[LucMeekes · Accepted solution]

With the rnd() random number generator you get random numbers from a uniform sitribution. The chance of hitting a specific number (such as 0) is virtually zero; exactly it is 1 over the number of possible random values, since Mathcad's range of real number is pretty large, you chance of hitting 0 therefor is pretty low: approximatly 0. If you want a specific probability of occurrence of any number, you have to restrict the number of possible random values. I admit that I restricted it fiercely. You may want to use:

Z(z,z)=: floor(rnd(n^2))/n^3

or any higher power

I considered another option to include 0, which is

Z(z,z)=: ((rnd(2*n)-n)^2)/n^2

This brings 0 to the middle of the range, but does not increase the probability of hitting it exactly. On first run with n=4, I got 3.33e-6 as smallest value.

Yet another possibility is to additionally declare all values within a certain distance from 0 as 0 itself:

Z(z,z)=:

              |r<-rnd(n)/n^2

              |return if(r<0.0001,0,r)

But this ruins your uniform distribution.

Another construct could be:

Z(z,z)=: if(rnd(n)<0.0001,0,rnd(r))/n^2

Which could result in non-zero numbers below 0.0001/n^2, because if on first call to rnd() the result is >= 0.0001, a new random number is generated by the second call to rnd(), but I don't exactly know what it does to your (uniform) distribution.

Success!
Luc

[anthonyQueen]

Thanks Luc again. I certainly try them too!!