mf_lusolve.hlp

是一个经济学管理应用软件 很难找的 但是经济学学生又必须用到

{smcl}
{* 28mar2005}{...}
{cmd:help mata lusolve()}
{hline}
{* index solve AX=B}{...}
{* index lusolve()}{...}
{* index _lusolve()}{...}
{* index _lusolve_la()}{...}
{* index LAPACK}{...}

{title:Title}

{p 4 8 2}
{bf:[M-5] lusolve() -- Solve AX=B for X using LU decomposition}


{title:Syntax}

{p 4 8 2}
{it:numeric matrix}
{cmd:lusolve(}{it:numeric matrix A}{cmd:,}
{it:numeric matrix B}{cmd:)}

{p 4 27 2}
{it:numeric matrix}
{cmd:lusolve(}{it:numeric matrix A}{cmd:,}
{it:numeric matrix B}{cmd:,} {it:real scalar tol}{cmd:)}


{p 4 8 2}
{it:void}{bind:         }
{cmd:_lusolve(}{it:numeric matrix A}{cmd:,}
{it:numeric matrix B}{cmd:)}

{p 4 27 2}
{it:void}{bind:         }
{cmd:_lusolve(}{it:numeric matrix A}{cmd:,}
{it:numeric matrix B}{cmd:,} {it:real scalar tol}{cmd:)}


{p 4 8 2}
{it:real scalar}{bind:  }
{cmd:_lusolve_la(}{it:numeric matrix A}{cmd:,}
{it:numeric matrix B}{cmd:)}

{p 4 8}
{it:real scalar}{bind:  }
{cmd:_lusolve_la(}{it:numeric matrix A}{cmd:,}
{it:numeric matrix B}{cmd:,} {it: real scalar tol}{cmd:)}



{title:Description}

{p 4 4 2}
{cmd:lusolve(}{it:A}{cmd:,} {it:B}{cmd:)} 
solves {it:A}{it:X}={it:B} and
returns {it:X}.  {cmd:lusolve()} returns a matrix of missing values if {it:A}
is singular.  

{p 4 4 2}
{cmd:lusolve(}{it:A}{cmd:,} {it:B}{cmd:,} {it:tol}{cmd:)}
does the same thing, but allows you to specify the tolerance for declaring 
that {it:A} is singular; see {it:Tolerance} under {it:Remarks} below.

{p 4 4 2}
{cmd:_lusolve(}{it:A}{cmd:,} {it:B}{cmd:)} 
and 
{cmd:_lusolve(}{it:A}{cmd:,} {it:B}{cmd:,} {it:tol}{cmd:)}
do the same thing except that, rather than returning the solution {it:X}, 
they overwrite {it:B} with the solution and, in the process of making the 
calculation, they destroy the contents of {it:A}.

{p 4 4 2}
{cmd:_lusolve_la(}{it:A}{cmd:,} {it:B}{cmd:)}
and 
{cmd:_lusolve_la(}{it:A}{cmd:,} {it:B}{cmd:,} {it:tol}{cmd:)} 
are the interfaces to the {bf:{help m1_lapack:[M-1] LAPACK}} routines that do
the work.  They solve {it:A}{it:X}=B for {it:X}, returning the solution in
{it:B} and, in the process, using as workspace (overwriting) {it:A}.
The routines return 1 if {it:A} was singular and 0 otherwise.  If {it:A} was
singular, {it:B} is overwritten with a matrix of missing values.


{title:Remarks}

{p 4 4 2}
The above functions solve {it:AX}={it:B} via LU decomposition and are very
accurate.  An alternative is 
{bf:{help mf_qrsolve:[M-5] qrsolve()}},
which uses QR decomposition.  The difference between
the two solutions is not, practically speaking, accuracy.  When {it:A} is of
full rank, both routines return equivalent results, and the LU approach is
quicker, using approximately {it:O}(2/3{it:n}^3) operations rather than
{it:O}(4/3{it:n}^3), where {it:A} is {it:n} {it:x} {it:n}.

{p 4 4 2}
The difference arises when {it:A} is singular.  In that case, the LU-based
routines documented here return missing values.  The QR-based routines
documented in {bf:{help mf_qrsolve:[M-5] qrsolve()}} return a generalized
(least-squares) solution.

{p 4 4 2}
For more information on LU and QR decomposition, see 
{bf:{help mf_lud:[M-5] lud()}}
and 
see {bf:{help mf_qrd:[M-5] qrd()}}.

{p 4 4 2}
Remarks are presented under the headings

	{bf:Derivation}
	{bf:Relationship to inversion}
	{bf:Tolerance}


{title:Derivation}

{p 4 4 2}
We wish to solve for {it:X}

		{it:A}{it:X} = {it:B}{right:(1)    }

{p 4 4 2}
Perform LU decomposition on {it:A} so that we have {it:A} =
{it:P}{it:L}{it:U}.  Then (1) can be written

		{it:P}{it:L}{it:U}{it:X} = {it:B}

{p 4 4 2}
or, premultiplying by {it:P}{bf:'} and remembering that
{it:P}{bf:'}{it:P}={it:I},

		{it:L}{it:U}{it:X} = P'{it:B}{right:(2)    }

{p 4 4 2}
Define 

		{it:Z} = {it:U}{it:X}{right:(3)    }

{p 4 4 2}
Then (2) can be rewritten

		{it:L}{it:Z} = {it:P}{bf:'}{it:B}{right:(4)    }

{p 4 4 2}
It is easy to solve (4) for {it:Z} because {it:L} is a lower-triangular
matrix.  Once {it:Z} is known, it is easy to solve (3) for {it:X} because
{it:U} is upper-triangular.



{title:Relationship to inversion}

{p 4 4 2}
Another way to solve

		{it:AX} = {it:B}

{p 4 4 2}
is to obtain {it:A}^(-1) and then calculate

		{it:X} = {it:A}^(-1){it:B}

{p 4 4 2}
It is, however, better to solve {it:AX}={it:B} directly because fewer
numerical operations are required, and the result is therefore more
accurate and obtained in less computer time.

{p 4 4 2}
Indeed, rather than thinking about how solve can be implemented via inversion,
it is more productive to think about how inversion can be implemented via
solve.  Obtaining {it:A}^(-1) amounts to solving

		{it:A}{it:X} = {it:I}

{p 4 4 2}
Thus {cmd:lusolve()} (or any other solve routine) can be used to 
obtain inverses.  The inverse of {it:A} can be obtained by coding 

		: {cmd:Ainv = lusolve(}{it:A}{cmd:, I(rows(}{it:A}{cmd:)))}

{p 4 4 2}
In fact, we provide {bf:{help mf_luinv:[M-5] luinv()}} for obtaining inverses
via LU decomposition, but {cmd:luinv()} amounts to making the above
calculation, although a little memory is saved because the matrix 
{it:I} is never constructed.

{p 4 4 2}
Hence, everything said about {cmd:lusolve()} applies equally to
{cmd:luinv()}.


{title:Tolerance}

{p 4 4 2}
The default tolerance used is 

		{it:eta} = (1e-13)*trace(abs({it:U}))/{it:n}

{p 4 4 2}
where {it:U} is the upper-triangular matrix of the LU decomposition of 
{it:A}: {it:n} {it:x} {it:n}.  {it:A} is declared to be singular if any
diagonal element of {it:U} is less than or equal to {it:eta}.

{p 4 4 2}
If you specify {it:tol}>0, the value you specify is used to multiply 
{it:eta}.  You may instead specify {it:tol}<=0, and then the negative of the
value you specify is used in place of {it:eta}; see 
{bf:{help m1_tolerance:[M-1] tolerance}}.

{p 4 4 2}
So why not specify {it:tol}=0?  
You do not want to do that because, as matrices become very close to being 
singular, results can become inaccurate.  Here is an example:

	: {cmd:uniformseed(12345)}

	: {cmd:A = lowertriangle(uniform(4,4))}

	: {cmd:A[3,3] = 1e-15}

	: {cmd:trux = uniform(4,1)}

	: {cmd:b = A*trux}

	: {cmd:/*} {it:the above created an Ax=b problem, and we have placed the true}
	:    {it:value of x in trux.  We now obtain the solution via lusolve()}
	:    {it:and compare trux with the value obtained:}
	: {cmd:*/}

	: {cmd:x = lusolve(A, b, 0)}

	: {cmd:trux, x}
	{res}       {txt}          1             2
	    {c TLC}{hline 29}{c TRC}
	  1 {c |}  {res}.7997150919   .7997150919{txt}  {c |}{...}
   <- {it:The discussed numerical}
	  2 {c |}  {res}.9102488109   .9102488109{txt}  {c |}{...}
      {it:instability can cause this}
	  3 {c |}  {res} .442547889   .3593109488{txt}  {c |}{...}
      {it:output to vary a little}
	  4 {c |}  {res} .756650276   .8337468202{txt}  {c |}{...}
      {it:across different computers}
	    {c BLC}{hline 29}{c BRC}

{p 4 4 2}
We would like to see the second column being nearly equal to the first -- 
the estimated {it:x} being nearly equal to the true {it:x} -- but there 
are substantial differences.

{p 4 4 2}
It is worth noting that, 
even through the difference between {cmd:x} and {cmd:xtrue} is substantial, 
the difference between them is small
in the prediction space:

	: {cmd:A*trux-b, A*x-b}
	{res}       {txt}           1              2
	    {c TLC}{hline 31}{c TRC}
	  1 {c |}  {res}           0   -2.77556e-17{txt}  {c |}
	  2 {c |}  {res}           0              0{txt}  {c |}
	  3 {c |}  {res}           0              0{txt}  {c |}
	  4 {c |}  {res}           0              0{txt}  {c |}
	    {c BLC}{hline 31}{c BRC}

{p 4 4 2}
What made this problem so difficult was the line {cmd:A[3,3] = 1e-15}.  Remove
that and you would find that the maximum difference between {cmd:x} and
{cmd:trux} would be 2.22045e-16.

{p 4 4 2}
The degree to which the residuals {bf:A*x-b} are a reliable measure of the
accuracy of {it:x} depends on the condition number of the matrix, which can be
obtained by {bf:{help mf_cond:[M-5] cond()}}, which in the case of {cmd:A}, is
3.2947e+15.  If the matrix is well conditioned (the condition number is
small), small residuals imply an accurate solution for {it:x}.  
If the matrix is ill conditioned, 
small residuals are not a reliable indicator of accuracy.

{p 4 4 2}
Another way to check the accuracy of {cmd:x} is to set {it:tol}=0 and to see
how well {it:x} could be obtained were {it:x}={cmd:x}:

	: {cmd:x  = lusolve(A, b, 0)}

	: {cmd:x2 = lusolve(A, A*x, 0)}

{p 4 4 2}
If {cmd:x} and {cmd:x2} are virtually the same, then you can safely assume 
that {cmd:x} is the result of a numerically accurate calculation.
You might compare {cmd:x} and {cmd:x2} using {cmd:mreldif(x2,x)}; see 
{bf:{help mf_reldif:[M-5] reldif()}}.  
In our example, {cmd:mreldif(x2,x)} is .03, a large difference.

{p 4 4 2}
The point of all of this is that if {it:A} is ill conditioned, then small
changes in {it:A} or {it:B} can lead to radical differences in the solutions
for {it:X}.


{title:Conformability}

    {cmd:lusolve(}{it:A}{cmd:,} {it:B}{cmd:,} {it:tol}{cmd:)}:
	{it:input:}
		{it:A}:  {it:n x n}
		{it:B}:  {it:n x k}
	      {it:tol}:  1 {it:x} 1    (optional)
	{it:output:}
	   {it:result}:  {it:n x k}

    {cmd:_lusolve(}{it:A}{cmd:,} {it:B}{cmd:,} {it:tol}{cmd:)}:
	{it:input:}
		{it:A}:  {it:n x n}
		{it:B}:  {it:n x k}
	      {it:tol}:  1 {it:x} 1    (optional)
	{it:output:}
		{it:A}:  0 {it:x} 0
		{it:B}:  {it:n x k}

    {cmd:_lusolve_la(}{it:A}{cmd:,} {it:B}{cmd:,} {it:tol}{cmd:)}:
	{it:input:}
		{it:A}:  {it:n x n}
		{it:B}:  {it:n x k}
	      {it:tol}:  1 {it:x} 1    (optional)
	{it:output:}
		{it:A}:  0 {it:x} 0
		{it:B}:  {it:n x k}
	   {it:result}:  1 {it:x} 1


{title:Diagnostics}

{p 4 4 2}
{cmd:lusolve(}{it:A}{cmd:,} {it:B}{cmd:,} ...{cmd:)},
{cmd:_lusolve(}{it:A}{cmd:,} {it:B}{cmd:,} ...{cmd:)},
and
{cmd:_lusolve(}_la{it:A}{cmd:,} {it:B}{cmd:,} ...{cmd:)}
return a result containing missing if 
{it:A} or {it:B} contain missing values.
The functions return a result containing all missing values if {it:A} is
singular.

{p 4 4 2}
{cmd:_lusolve(}{it:A}{cmd:,} {it:B}{cmd:,} ...{cmd:)}
and 
{cmd:_lusolve_la(}{it:A}{cmd:,} {it:B}{cmd:,} ...{cmd:)}
abort with error if {it:A} or {it:B} is a view.

{p 4 4 2}
{cmd:_lusolve_la(}{it:A}{cmd:,} {it:B}{cmd:,} ...{cmd:)}
should not be used directly; use 
{cmd:_lusolve()}.


{title:Source code}

{p 4 4 2}
{view lusolve.mata, adopath asis:lusolve.mata},
{view _lusolve.mata, adopath asis:_lusolve.mata};
{cmd:_lusolve_la()} is built-in.


{title:Also see}

{p 4 13 2}
Manual:  {hi:[M-5] lusolve()}

{p 4 13 2}
Online:  help for 
{bf:{help mf_luinv:[M-5] luinv()}},
{bf:{help mf_lud:[M-5] lud()}},
{bf:{help mf_solvelower:[M-5] solvelower()}},
{bf:{help mf_cholsolve:[M-5] cholsolve()}},
{bf:{help mf_qrsolve:[M-5] qrsolve()}},
{bf:{help mf_svsolve:[M-5] svsolve()}};
{bf:{help m4_matrix:[M-4] matrix}}
{p_end}