Inverse of matrix

[Clara Etcheverry]

Hi !
I have to find the inverse of a matrix 262*262
When I use qrinv and pinv it didn't give the same result, is it normal?
Moreover, the function inv didn't work...
Clara

[Jean-Claude Arbaut]

There is no inv function in Mata, but there is luinv, qrinv and pinv, and also invsym for symmetric pos. def. matrices. As the names suggest, luinv uses the LU decomposition, qrinv uses the QR decomposition and pinv is the pseudoinverse and uses the SVD decomposition. The last, invsym, probably uses the Cholesky decomposition, which is similar to the LU decomposition.

The computations are thus different, and even if mathematically they are supposed to yield the same result, the actual numerical result will vary due to the finite precision of floating-point computations. In good cases, the norm of the difference will be small, but numerical instability with badly conditioned matrices can yield very different inverses. In such cases, I would use the pseudoinverse, which should be more robust (LU is the least robust, QR is usually a good compromise).

Given that the machine epsilon is 2.22e-16 with the usual "double" precision, and that the norm of the difference is going to increase with the dimension of the matrix, a value of the order 1e-12 to 1e-10 should not be considered a problem here.
If the difference is much larger, you matrix may be badly conditioned (common with multicollinearity), and you should investigate further: that could me removing variables, switching to a method more robust to multicollinearity...


Here is an example:

set seed 17760704
mata
a = rnormal(262,262,0,1)
u = luinv(a)
v = qrinv(a)
w = pinv(a)
norm(u-v)
norm(v-w)
end

For this case I get respectively 7.25174e-11 and 7.48100e-11.

For more informations about finite-precision arithmetic, I suggest this excellent article by Samuel Goldberg:
What Every Computer Scientist Should Know About Floating-Point Arithmetic
http://www.validlab.com/goldberg/paper.pdf

Hope this helps

Jean-Claude Arbaut

[Clara Etcheverry]

Thank you for your answer!!
The different is much larger with my matrix... Does it exist more robust method?
Moreover, with pinv and qring the matrix are very very different!
Clara

[Clara Etcheverry]

(The order is weird, it is equal to 410.4698504)

[Jean-Claude Arbaut]

What's the value of the relative error, norm(a-b)/norm(a)? (with a and b the matrix inverses)

And where do this matrix come from? Could you describe a little bit how it's built?

[Clara Etcheverry]

norm (a-b)/norm a = 0,0446224955271596

The matrix is (I*S*I)^-1 where I is the identity matrix262*262 and S is a matrix 262*262

Thank you very much!

[Jean-Claude Arbaut]

If I is the identity, I*S*I = S, so your expression simplifies to S^-1. What is S here?

[Clara Etcheverry]

I am modelling the demand and the supply for a market. I have 262 alternatives that the consumer can choose (that's why there are matrix of dimension 262*262).
And sorry I forgot to precise, I is the (262x262) ownership matrix with element 1 if the product j is sold by the firm f and 0 otherwise.

The matrix S is a matrix 262*262 where

Where (but really long and difficult to understand....)

egen pair=group(RetailGroup Brand)
sum pair
global nb_pair=r(max)
mkmat pair if mois==1

mata somme2=J($nb_p*$T,$nb_simul,1)
forvalues p=1/$nb_p {
mata somme2`p'=somme
forvalues i=1/$nb_p {
if pair[`p',1]==pair[`i',1] {
mata somme2`p'=somme2`p'-totob`i'
}
}
}


forvalues k=1/$nb_p {
forvalues j=1/$nb_p {
if pair[`k',1]!=pair[`j',1] {
mata temp`k'_`j'=totob`k':/somme2`j'- totob`k':/somme
}
else if pair[`k',1]==pair[`j',1] {
mata temp`k'_`j'=totob`k':/somme
}
}
}

/*market shares*/

forvalues k=1/$nb_p {
mata : mata drop totob`k'
mata S`k'=J($nb_p*$T,1,0)
forvalues j=1/$nb_p {
mata part_sim`k'_`j'=(temp`k'_`j'*J($nb_simul,1,(1/$nb_simul)))
mata : mata drop temp`k'_`j'
sort mois produit

/* k vecteurs with market share and differences of market share/*
mata S`k'=S`k'+part_sim`k'_`j':*p`j'
mata : mata drop part_sim`k'_`j'
}
}

/*I stick the vectors*/
mata S=J($nb_p*$T,1,0)

forvalues k=1/$nb_p {
mata S=S,S`k'
mata : mata drop S`k'
}
mata st_matrix("S", S)

/* l'enleve la colonne avec que des 0*/
mata S=S[.,2..$nb_p+1]


/*S for each month

forvalues t=1/$T {
mata ST`t'=J($nb_p,$nb_p,0)

mata ST`t'=S[(`t'-1)*$nb_p+1..`t'*$nb_p,1...]

mata st_matrix("ST`t'", ST`t')
}



(nb_p number of alternative and T the period).


Clara