Inverting a large, dense matrix

In this post I show some experiences with inverting a large, dense matrix. The variability in performance was a bit surprising to me. 

Background

I have available in a GAMS GDX file a large matrix \(A\). This is an input-output matrix. The goal is to find the so-called Leontief inverse \[(I-A)^{-1}\] and return this in a GAMS GDX for subsequent use [1]. It is well-known that instead of calculating the inverse it is better to form an LU factorization of \(I-A\) and then use that for different rhs vectors \(b\) (essentially a pair of easy triangular solves). [2] However, economists are very much used to having access to an explicit inverse.

Expanding Visual Studio in the Task Manager

 

This is scary. Note that is in an idle state (nothing compiling, running, or debugging at the moment). Soon we all need 128 GB of RAM.

GAMS: Undocumented PUT formats

The GAMS documentation mentions f.nr=1 (standard notation) and f.nr=2 (scientific notation) for formatting numeric items using the PUT statement. There are however a few more. Here we print values using different formats:

       f.nr=1       f.nr=2       f.nr=3       f.nr=4
         1.20     1.20E+00 1.2000000000          1.2
         1.23     1.23E+00 1.2345000000       1.2345
 123450000.00     1.23E+08 123450000.00    123450000
         0.00     1.20E-07 0.0000001200       1.2E-7
         0.00     1.23E-07 0.0000001235    1.2345E-7

It looks like formats 3 and 4 are inspired by SAS PUT formats f12.10 and best12.10.

References

1. Dictionary of formats, https://documentation.sas.com/doc/en/pgmsascdc/9.4_3.5/leforinforref/p0z62k899n6a7wn1r5in6q5253v1.htm

Appendix: GAMS code

set    i   'values'   /i1*i5/    fmt 'formats'  /'f.nr=1'*'f.nr=4'/ ; parameter p(i) 'test values' /    i1  1.2    i2  1.2345    i3  1.2345e8    i4  1.2e-7    i5  1.2345e-7 /; * * write all values in all formats * file f /put.txt/; * right align text f.lj = 1; put f; * write headers loop(fmt,    put fmt.tl:13 ); put /; * write values loop(i,    loop(fmt,       put ' ':0;       f.nr=ord(fmt);       put p(i);    );    put /; );

Rewriting old GAMS code

It is always a good idea to revisit existing GAMS code and see if we can improve it. Here is an example of an actual model.

The problem is that we want to set up a mapping set between two sets based on the first two characters. If they are the same, the pair should be added to the mapping. The old code looked like:

sets    r 'regions' /        AP   'Appalachia'        CB   'Corn Belt'        LA   'Lake States'        /    s 'subregions' /        APN0501        APN0502        APN0504        CBL0704        CBM0404        LAM0404        LAM0406        /    rs(r,s) 'mapping regions<->subregions' ; * old code: PARAMETER value; LOOP(R,  value('1',R) = 100*ord( R.tl,1) + ord( R.tl,2); ); LOOP(S,  value('2',S) = 100*ord( S.tl,1) + ord( S.tl,2); ); RS(R,S)$(value('1',R) = value('2',S)) = YES; display value; display rs;

  

To verify this indeed generates the correct result for set rs, we look at the output:

----     27 PARAMETER value  

           AP          CB          LA     APN0501     APN0502     APN0504     CBL0704     CBM0404     LAM0404

1    6580.000    6766.000    7665.000
2                                        6580.000    6580.000    6580.000    6766.000    6766.000    7665.000

+     LAM0406

2    7665.000


----     28 SET rs  mapping regions<->subregions

       APN0501     APN0502     APN0504     CBL0704     CBM0404     LAM0404     LAM0406

AP         YES         YES         YES
CB                                             YES         YES
LA                                                                     YES         YES

The set rs is obviously correctly populated. The auxiliary parameter value contains the ASCII values: e.g. AP has value 6580 as the decimal ASCII codes for A and P are 65 and 80.

Notes:

• The function ord(string, pos) returns the ASCII value of the character at position pos of the string.
• The suffix .TL is the text string of the set element (internally GAMS does not work with the strings but rather an integer id for each set element; here we make it explicit we want the string).
• The loops are not really needed. We could have used:
  value('1',R) = 100*ord( R.tl,1) + ord( R.tl,2);
  value('2',S) = 100*ord( S.tl,1) + ord( S.tl,2);
• The mapping set is a very powerful concept in GAMS. It is somewhat like a dict in Python, except it is not one way. It can map from r to s but also from s to r.

We can simplify the code to just:

* new code: * compare first two characters in r and s rs(r,s) = ord(r.tl,1)=ord(s.tl,1) and ord(r.tl,2)=ord(s.tl,2); display rs;

This gives the same result:

----     22 SET rs  mapping regions<->subregions

       APN0501     APN0502     APN0504     CBL0704     CBM0404     LAM0404     LAM0406

AP         YES         YES         YES
CB                                             YES         YES
LA                                                                     YES         YES

Notes:

• No need for the auxiliary parameter value.
• No loops.
• We could add a check that we have no unmapped subregions. E.g. by:
  abort$(card(rs)<>card(s)) "There are unmapped subregions";
• This will also catch cases where the casing does not match.
• Critical note: GAMS has very poor support for strings. This should have been rectified a long time ago. There is no excuse for this deficiency. 

Conclusion

It is always worthwhile to revisit existing models and clean them up. During development, especially when under time pressure, we often are happy once our code works. As the lifetime of models can be surprisingly long, and thus new generations of users and modelers start working with old models, it is a good idea to do a cleanup round to make the code as readable as possible.