0

How can I generate a random triangular matrix? (upper and lower)

Usually I use rand(n) but if I try tril(rand(n)) it will be singular and I don't want that.

Giuseppe Trapasso
  • 189
  • 2
  • 14

3 Answers3

3

your answer is correct:

 A=tril(rand(n))

you can check that this matrix is not singular using 

 rcond(A)>eps

or

 min(svd(A))>eps

and verifying that the smallest singular value is larger than eps, or any other numerical tolerance that is relevant to your needs. (the code will return 1 or 0). For n>50 you'll start to approach singular matrices.

Here's a small analysis of how the matrix approaches singularity with it's size...

enter image description here

bla
  • 25,846
  • 10
  • 70
  • 101
  • 1
    Wow! Faster than exponential behavior. Wondering why? – Severin Pappadeux Oct 07 '15 at 17:34
  • Actually I didn'the want a ill-conditionated matrix either, ando my 'n' was 100. That could be the problem... – Giuseppe Trapasso Oct 07 '15 at 20:48
  • Do similar matrices have the same Reciprocal condition number? – BillBokeey Oct 08 '15 at 12:34
  • Triangular matrices are nonsingular if there is no zero on the diagonal by definition. You can get the same problems with `tril(rand(100)+100)` this is a genuine implementation problem and this analysis cannot be trusted. – percusse Oct 08 '15 at 12:59
  • @percusse `matrices are nonsingular if there is no zero on the diagonal by definition` Nonsense! Matrix is nonsingular if determinant is non-zero. I could make singular tri-matrix with non-zero diagonal with easy. – Severin Pappadeux Oct 08 '15 at 15:08
  • @severin Please do one for the fun of it. I'll be waiting. – percusse Oct 08 '15 at 15:24
  • @SeverinPappadeux And anything yet? – percusse Oct 09 '15 at 05:20
  • @SeverinPappadeux Still nothing? And until you prove that nonsense comment. – percusse Oct 10 '15 at 01:07
  • There is a big difference between matrices that are close from a singularity with respect to a certain norm and matrices that have no singularity in their neighborhood.. It's the definition of conditioning. I of course agree with you percusse, that a matrix will be singular if and only if it has a zero on its diagonal, but you can still be very close from a singularity (very ill-conditionated then), and have a non zero determinant. That's the point in @bla answer – BillBokeey Oct 12 '15 at 11:37
  • @BillBokeey Not with triangular matrices. Clearly `tril(rand(100)+0.1)` is nonsingular. But not according to its condition number (run a few times). It is a geniune analysis error which is related how determinant and eigenvalues are computed for such matrices. At those sizes, you need to preprocess your matrices first or call the proper Fortran routine that is meant for triangular matrices to assess numerical stability. You cannot just go `det(A)` if you don't know what your matrix looks like which is the problem here. – percusse Oct 12 '15 at 14:12
  • @BillBokeey The example I gave is not close to any singular matrix. But if you check its condition number you will get the impression that it is. So no we are not saying the same thing. The analysis in that plot is wrong. It shows the success of the rcond and svd analysis which gets worse and worse without matrix getting singular at all. – percusse Oct 12 '15 at 20:47
0

Ok, here are the thoughts on the triangular matrix singularity. Determinant of the triangular matrix is what determines singularity, because it is going into denominator when inverse matrix is build. The property of the triangular matrix is such that determinant equals to the product of diagonal elements.

So for matrix NxN we on the diagonal have the product of i.i.d. U(0,1) numbers. Obviously, determinant will decrease as N increases due to the fact that all number are <1 and more you have, the less value the product (aka determinant) will have.

It is interesting to check that for det=X1X2...*XN the mean value will go down as 2-N, because each term in product is U(0,1) with mean of 1/2, and they are all i.i.d. Alternative check would be to compute mean from product PDF (see https://math.stackexchange.com/questions/659254/product-distribution-of-two-uniform-distribution-what-about-3-or-more), and, indeed, it will give you exactly the same result, 2-N. Variance of the determinant could be computed as well, as second momentum minus mean squared, and it is equal to (3-N-4-N).

Please note, those are mean values, which you could expect on average, say, if you sample 106 triangular matrices with N=100, compute their determinant and average it, you should find it to be pretty close to 2-100.

That's where the problem lies. On average, triangular random matrix goes to singularity exponentially with the grows of N. 2-10 is about equal to 1/1,000. 2-20 is about equal to 1/1,000,000. For N=100 it should be, on average, around 10-30 or so, which makes whole exercise moot.

Unfortunately, I cannot offer anything beyond this simple analysis.

Community
  • 1
  • 1
Severin Pappadeux
  • 18,636
  • 3
  • 38
  • 64
0

If you want a well-conditioned random triangular matrix you could take the triangular part of A2 with A = rand(n). So triu(A * A) for any size n is well conditioned but of course has complexity O(n3) for the matrix-matrix multiplication.

ps3lister
  • 47
  • 9