Questions about lud() and lusolve()

[Anna Polsi]

Hi there!

I have several questions about lud()and lusolve() functions.

Reading MATA manual relative to LU factorization - i.e., lud(numeric matrix A, L, U, p) function - I do not really get the arguments Mata is expected to receive other than the known matrix A because L and U come from LU decomposition - so, they are not yet known. Also, I do not understand how the permutation vector p is chosen. What confused me even more is the statement in the manual "2. The types of L, U, p, and q are irrelevant; results are returned there." because if I do not create L, U, and p I get an error because Mata does not find them, but if I create them with, for example, zero or empty elements, they are not filled after the decomposition.
I also tried to use the function _lud_la(numeric matrix A, q) setting q = (1\2\2) as documented in MATA manual. The function performs the LU decomposition smoothly but the original matrix A is modified from which L and U matrices can be extrapolated as explained in the documentation. But I do not want the original matrix A to be modified that is why I want to use lud(). So, how am I supposed to proceed?

Finally, I have a computational doubt: Is it correct that cholsolve() and lusolve() return the same solution (i.e., a vector x) when the matrix A is symmetric (and, of course, positive definite and square)? For example,

A = (2,-1\-1,2) 
b = J(rows(A),1,1)
MatrixLUSolve = lusolve(A,b) 
MatrixCholSolve = cholsolve(A, b)

Here, the solution is the same even though AX=B is solved for X using two different decompositions.



Thank you,

Anna

[William Lisowski]

This example demonstrates that the names L, U, and p can be assigned as real scalars with an arbitrary value, and lud() returns the results to those names, ignoring their original type.

. mata:
------------------------------------------------- mata (type end to exit) ----------------------
: rseed(12345)

: A = runiform(4,4)

: L = 0

: U = 0

: p = 0

: mata describe

      # bytes   type                        name and extent
-------------------------------------------------------------------------------
          128   real matrix                 A[4,4]
            8   real scalar                 L
            8   real scalar                 U
            8   real scalar                 p
-------------------------------------------------------------------------------

: lud(A, L, U, p)

: mata describe

      # bytes   type                        name and extent
-------------------------------------------------------------------------------
          128   real matrix                 A[4,4]
          128   real matrix                 L[4,4]
          128   real matrix                 U[4,4]
           32   real colvector              p[4]
-------------------------------------------------------------------------------

: A
                 1             2             3             4
    +---------------------------------------------------------+
  1 |  .3576297229    .400442617    .689383317   .5597355706  |
  2 |   .574451294   .2076905269   .0286626993   .6889244779  |
  3 |  .4693433591   .2071525958   .0039322524   .0130297115  |
  4 |  .4204224235   .6161914501    .894883612   .4106360552  |
    +---------------------------------------------------------+

: L
                  1              2              3              4
    +-------------------------------------------------------------+
  1 |             1              0              0              0  |
  2 |   .7318678327              1              0              0  |
  3 |   .6225588255    .5841215386              1              0  |
  4 |   .8170289875    .0807071685   -.5588601841              1  |
    +-------------------------------------------------------------+

: U
                  1              2              3              4
    +-------------------------------------------------------------+
  1 |    .574451294    .2076905269    .0286626993    .6889244779  |
  2 |             0    .4641894344    .8739063044   -.0935656093  |
  3 |             0              0    .1610716055    .1854932445  |
  4 |             0              0              0    -.438625353  |
    +-------------------------------------------------------------+

: p
       1
    +-----+
  1 |  3  |
  2 |  1  |
  3 |  4  |
  4 |  2  |
    +-----+

: end
------------------------------------------------------------------------------------------------

Is it correct that cholsolve() and lusolve() return the same solution (i.e., a vector x) when the matrix A is symmetric (and, of course, positive definite and square)?

Both functions solve the same equation AX = B, so when they both can solve the equation, they both must return the same solution.

[Anna Polsi]

Thank you for the clarification.