Skip to main content
A
anthonyQueen11-Garnet
11-Garnet
July 11, 2019
Solved

How to generate a semipositive indecomposible matrix of nxn

  • July 11, 2019
  • 1 reply
  • 2131 views

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.

 

Best answer by LucMeekes

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

          

1 reply

L
LucMeekes23-Emerald IV
23-Emerald IV
July 11, 2019
Use Z (z,z)=:floor (rnd (n))/n^2.
That will cause a 0 every now and then.

Success!
Luc
    A
    anthonyQueen11-GarnetAuthor
    11-Garnet
    July 11, 2019

    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. 

    L
    LucMeekes23-Emerald IVAnswer
    23-Emerald IV
    July 12, 2019

    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

              

      Terms & ConditionsCookie settingsAccessibility statement