mf_qrd.hlp

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

{smcl}
{* 31mar2005}{...}
{cmd:help mata qrd()}
{hline}
{* index QR decomposition}{...}
{* index qrd()}{...}
{* index hqrd()}{...}
{* index _hqrd()}{...}
{* index hqrdmultq()}{...}
{* index hqrdmultq1t()}{...}
{* index hqrdr()}{...}
{* index hqrdr1()}{...}
{* index qrdp()}{...}
{* index hqrdp()}{...}
{* index _hqrdp()}{...}
{* index _hqrdp_la()}{...}
{* index LAPACK}{...}
{* index decomposition}{...}

{title:Title}

{p 4 8 2}
{bf:[M-5] qrd() -- QR decomposition}


{title:Syntax}

{p 8 12 2}
{it:void}{bind:           }
{cmd:qrd(}{it:numeric matrix A}{cmd:,} {it:Q}{cmd:,} {it:R}{cmd:)}


{p 8 12 2}
{it:void} {bind:          }
{cmd:hqrd(}{it:numeric matrix A}{cmd:,} {it:H}{cmd:,} {it:tau}{cmd:,} {it:R1}{cmd:)}

{p 8 12 2}
{it:void}{bind:          }
{cmd:_hqrd(}{it:numeric matrix A}{cmd:,} {it:tau}{cmd:,} {it:R1}{cmd:)}



{p 8 34 2}
{it:numeric matrix}{bind: }
{cmd:hqrdmultq(}{it:numeric matrix H}{cmd:,} {it:numeric matrix tau}{cmd:,}
	{it:numeric matrix X}{cmd:,}  {it:real scalar transpose}{cmd:)}

{p 8 36 2}
{it:numeric matrix}{bind: }
{cmd:hqrdmultq1t(}{it:numeric matrix H}{cmd:,} {it:numeric matrix tau}{cmd:,}
	{it:numeric matrix X}{cmd:)}

{p 8 36 2}
{it:numeric matrix}{bind: }
{cmd:hqrdq(}{it:numeric matrix H}{cmd:,} {it:numeric matrix tau}{cmd:)}

{p 8 12 2}
{it:numeric matrix}{bind: }
{cmd:hqrdq1(}{it:numeric matrix H}{cmd:,} {it:numeric matrix tau}{cmd:)}

{p 8 12 2}
{it:numeric matrix}{bind: }
{cmd:hqrdr(}{it:numeric matrix H}{cmd:)}

{p 8 12 2}
{it:numeric matrix}{bind: }
{cmd:hqrdr1(}{it:numeric matrix H}{cmd:)}



{p 8 12 2}
{it:void}{bind:           }
{cmd:qrdp(}{it:numeric matrix A}{cmd:,} {it:Q}{cmd:,} {it:R}{cmd:,} 
	{it:real rowvector p}{cmd:)}

{p 8 12 2}
{it:void}{bind:           }
{cmd:hqrdp(}{it:numeric matrix A}{cmd:,} {it:H}{cmd:,} {it:tau}{cmd:,} {it:R1}{cmd:,}
	{it:real rowvector p}{cmd:)}

{p 8 12 2}
{it:void}{bind:          }
{cmd:_hqrdp(}{it:numeric matrix A}{cmd:,} {it:tau}{cmd:,} {it:R1}{cmd:,}
	{it:real rowvector p}{cmd:)}


{p 8 12 2}
{it:void}{bind:          }
{cmd:_hqrdp_la(}{it:numeric matrix A}{cmd:,} {it:tau}{cmd:,} {it:real rowvector p}{cmd:)}


{title:Description}

{p 4 4 2}
{cmd:qrd(}{it:A}{cmd:,} {it:Q}{cmd:,} {it:R}{cmd:)} 
calculates the QR decomposition of {it:A}: {it:m} {it:x} {it:n},
{it:m}>={it:n}, returning results in {it:Q} and {it:R}.

{p 4 4 2}
{cmd:hqrd(}{it:A}{cmd:,} {it:H}{cmd:,} {it:tau}{cmd:,} {it:R1}{cmd:)}
calculates the QR decomposition of {it:A}: {it:m} {it:x} {it:n},
{it:m}>={it:n}, but rather than returning {it:Q} and {it:R}, returns the
Householder vectors in {it:H} and the scale factors {it:tau} -- from which
{it:Q} can be formed -- and returns an upper-triangular matrix in {it:R1} that
is a submatrix of {it:R}; see {bf:Remarks} below for its definition.  Doing
this saves calculation and memory, and other routines allow you to manipulate
these matrices:

{p 8 12 2}
1.
    {cmd:hqrdmultq(}{it:H}{cmd:,} {it:tau}{cmd:,} {it:X}{cmd:,} {it:transpose}{cmd:)}
    returns {it:Q}{it:X} or {it:Q}{bf:'}{it:X} 
    based on the {it:Q} implied by {it:H} and {it:tau}.
    {it:Q}{it:X} is returned if
    {it:transpose}==0, and {it:Q}{bf:'}{it:X} is returned otherwise.

{p 8 12 2}
2.
    {cmd:hqrdmultq1t(}{it:H}{cmd:,} {it:tau}{cmd:,} {it:X}{cmd:)}
    returns {it:Q1}{bf:'}{it:X} 
    based on the {it:Q1} implied by {it:H} and {it:tau}.

{p 8 12 2}
3.
    {cmd:hqrdq(}{it:H}{cmd:,} {it:tau}{cmd:)} 
    returns the {it:Q} matrix implied by {it:H} and {it:tau}.  This function
    is rarely used.

{p 8 12 2}
4.
    {cmd:hqrdq1(}{it:H}{cmd:,} {it:tau}{cmd:)} 
    returns the {it:Q1} matrix implied by {it:H} and {it:tau}.  This function
    is rarely used.

{p 8 12 2}
5.
    {cmd:hqrdr(}{it:H}{cmd:)} 
    returns the full {it:R} matrix.  This function is rarely used.
    (It may surprise you that {cmd:hqrdr()} is a function of {it:H} and 
    not {it:R1}.  {it:R1} also happens to be stored in {it:H}, and there 
    is other useful information there, as well.)

{p 8 12 2}
6.
    {cmd:hqrdr1(}{it:H}{cmd:)} 
    returns the {it:R1} matrix.  This function is rarely used.

{p 4 4 2}
{cmd:_hqrd(}{it:A}{cmd:,} {it:tau}{cmd:,} {it:R1}{cmd:)} 
does the same thing as 
{cmd:hqrd(}{it:A}{cmd:,} {it:H}{cmd:,} {it:tau}{cmd:,} {it:R1}{cmd:)},
except that it overwrites {it:H} into {it:A} and so conserves even more
memory.

{p 4 4 2}
{cmd:qrdp(}{it:A}{cmd:,} {it:Q}{cmd:,} {it:R}{cmd:,} {it:p}{cmd:)} 
is similar to 
{cmd:qrd(}{it:A}{cmd:,} {it:Q}{cmd:,} {it:R}):  
it returns the QR
decomposition of {it:A} in {it:Q} and {it:R}.  The difference is that this
routine allows for pivoting.  New argument {it:p} specifies whether a column
is available for pivoting and, on output, {it:p} is overwritten with a
permutation vector that records the pivoting actually performed.  On input,
{it:p} can be specified as {cmd:.} (missing) -- meaning all columns are
available for pivoting -- or {it:p} can be specified as an {it:n} {it:x} 1
column vector containing 0s and 1s, with 1 meaning the column is fixed and so
may not be pivoted.

{p 4 4 2}
{cmd:hqrdp(}{it:A}{cmd:,} {it:H}{cmd:,} {it:tau}{cmd:,} {it:R1}{cmd:,} {it:p}{cmd:)} 
is a generalization of 
{cmd:hqrd(}{it:A}{cmd:,} {it:H}{cmd:,} {it:tau}{cmd:,} {it:R1}{cmd:)} 
just as {cmd:qrdp()} is a generalization of {cmd:qrd()}.

{p 4 4 2}
{cmd:_hqrdp(}{it:A}{cmd:,} {it:tau}{cmd:,}  {it:R1}{cmd:,} {it: p}{cmd:)} 
does the same thing as 
{cmd:hqrdp(}{it:A}{cmd:,} {it:H}{cmd:,} {it:tau}{cmd:,} {it:R1},
{it:p}{cmd:)}, except that {cmd:_hqrdp()} overwrites {it:H} into {it:A}.

{p 4 4 2}
{cmd:_hqrdp_la()} is the interface into the 
{bf:{help m1_lapack:[M-1] LAPACK}} routine that performs the QR calculation;
it is used by all the above routines.  Direct use of {cmd:_hqrdp_la()} is not
recommended.


{title:Remarks}

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

	{bf:QR decomposition}
	{bf:Avoiding calculation of Q}
	{bf:Pivoting}
	{bf:Least-squares solutions with dropped columns}


{title:QR decomposition}

{p 4 4 2}
The decomposition of square or nonsquare matrix {it:A} can be written 

			{it:A} = {it:QR}{right:(1)   }

{p 4 4 2}
where {it:Q} is an orthogonal matrix ({it:Q}{cmd:'}{it:Q} = {it:I}), and
{it:R} is upper triangular.  
{cmd:qrd(}{it:A}{cmd:,} {it:Q}{cmd:,} {it:R}{cmd:)}
 will make this calculation:

	: {cmd:A}
	{res}       {txt}1   2
	    {c TLC}{hline 9}{c TRC}
	  1 {c |}  {res}7   4{txt}  {c |}
	  2 {c |}  {res}9   6{txt}  {c |}
	  3 {c |}  {res}9   6{txt}  {c |}
	  4 {c |}  {res}7   2{txt}  {c |}
	  5 {c |}  {res}3   1{txt}  {c |}
	    {c BLC}{hline 9}{c BRC}

	: {cmd:Q = R = .}

	: {cmd:qrd(A, Q, R)}

	: {cmd:Ahat = Q*R}

	: {cmd:mreldif(Ahat, A)}
	  3.55271e-16


{title:Avoiding calculation of Q}

{p 4 4 2}
In fact, you probably do not want to use {cmd:qrd()}.  Calculating the 
necessary ingredients for {it:Q} is not too difficult, but going from 
those necessary ingredients to form {it:Q} is devilish.  In most cases, 
the necessary ingredients are all you need.  Those necessary ingredients 
are the Householder vectors and their scale factors, known as {it:H} and
{it:tau}.  For instance, one can write down a mathematical function
{it:f}({it:H}, {it:tau}, {it:X}) that will calculate {it:Q}{it:X} or
{it:Q}{cmd:'}{it:X} for some matrix {it:X}.

{p 4 4 2}
In addition, QR decomposition is often carried out on violently nonsquare
matrices {it:A}: {it:m} {it:x} {it:n}, {it:m}>>{it:n}.  We can write


					       {c TLC}    {c TRC}
		           {c TLC}               {c TRC}   {c |} {it:R1} {c |}
		A       =  {c |}  {it:Q1}      {it:Q2}   {c |}   {c |}    {c |}  =  {it:Q1}*{it:R1}  + {it:Q2}*{it:R2}
	                   {c BLC}               {c BRC}   {c |} {it:R2} {c |}      
					       {c BLC}    {c BRC}
	      {it:m x n         m x n   m x m-n     n x n     m x n    m x n}
                                              {it:m-n x n}


{p 4 4 2}
It turns out that {it:R2} is zero, and thus

					       {c TLC}    {c TRC}
		           {c TLC}               {c TRC}   {c |} {it:R1} {c |}
		A       =  {c |}  {it:Q1}      {it:Q2}   {c |}   {c |}    {c |}  =  {it:Q1}*{it:R1}
	                   {c BLC}               {c BRC}   {c |}  0 {c |}      
					       {c BLC}    {c BRC}
	      {it:m x n         m x n   m x m-n     n x n     m x n}
                                              {it:m-n x n}


{p 4 4 2}
Thus it is enough to know {it:Q1} and {it:R1}.  Rather than defining QR
decomposition as

			{it:A} = {it:QR},            {it:Q}: {it:m x m},          {it:R}: {it:m x n}{right:(1)   }

{p 4 4 2}
it is better to define it as 

			{it:A} = {it:Q1}*{it:R1}         {it:Q1}: {it:m x n}          {it:R1}: {it:n x n}{right:(1')  }

{p 4 4 2}
To appreciate the savings, consider the reasonable case where {it:m}=4000 and
{it:n}=3:

			{it:A} = {it:QR},            {it:Q}: 4000 {it:x} 4000,     {it:R}: 4000 {it:x} 3
    versus, 
			{it:A} = {it:Q1}*{it:R1}         {it:Q1}: 4000 {it:x} 3        {it:R1}:    3 {it:x} 3

{p 4 4 2}
Memory consumption is reduced from 125,094 kilobytes to 94 kilobytes, a 99.92
percent saving!

{p 4 4 2}
Combining the arguments,
we need not save {it:Q} because {it:Q1} is sufficient, 
we need not calculate {it:Q1} because {it:H} and {it:tau} are sufficient, and
we need not store {it:R} because {it:R1} is sufficient.

{p 4 4 2}
That is what 
{cmd:hqrd(}{it:A}{cmd:,} {it:H}{cmd:,} {it:tau}{cmd:,} {it:R1}{cmd:)}
does.  Having
used {cmd:hqrd()}, if you need to multiply the full {it:Q} by some matrix
{it:X}, you can use {cmd:hqrdmultq()}.  Having used {cmd:hqrd()}, if you need
the full {it:Q}, you can use {cmd:hqrdq()} to obtain it, but by that point you
will be making the devilish calculation you sought to avoid and so you might
as well have used {cmd:qrd()} to begin with.  If you want {it:Q1}, you can use
{cmd:hqrdq1()}.  Finally, having used {cmd:hqrd()}, if you need {it:R} or
{it:R1}, you can use {cmd:hqrdr()} and {cmd:hqrdr1()}:

	: {cmd:A}
	{res}       {txt}1   2
	    {c TLC}{hline 9}{c TRC}
	  1 {c |}  {res}7   4{txt}  {c |}
	  2 {c |}  {res}9   6{txt}  {c |}
	  3 {c |}  {res}9   6{txt}  {c |}
	  4 {c |}  {res}7   2{txt}  {c |}
	  5 {c |}  {res}3   1{txt}  {c |}
	    {c BLC}{hline 9}{c BRC}

	: {cmd:H = tau = R1 = .}

	: {cmd:hqrd(A, H, tau, R1)}

	: {cmd:Ahat = hqrdq1(H, tau) * R1}{...}
{col 55}// i.e., {it:Q1}*{it:R1}

	: {cmd:mreldif(Ahat, A)}
	  3.55271e-16


{title:Pivoting}

{p 4 4 2}
The QR decomposition with column pivoting solves