Re: Matrix norms? ..
- To: mathgroup at smc.vnet.net
- Subject: [mg6793] Re: Matrix norms? ..
- From: Daniel Lichtblau <danl>
- Date: Sun, 20 Apr 1997 01:50:32 -0400 (EDT)
- Organization: Wolfram Research, Inc.
- Sender: owner-wri-mathgroup at wolfram.com
Chung Cheong Ming wrote:
>
> Hello everyone,
>
> I just wonder if Mathematica can find the p-norms of a matrix?
> p can be 1, 2, infinity, though I know I can program those norms.
>
> Joseph
>
The infinity-norm is just the maximum over all rows of the sum of
absolute values. The 1-norm is the same, summed over columns.
Alternatively, it can be viewed as the infinity-norm of the transposed
matrix. The 2-norm is just the largest singular value.
The examples below were run in version 3 of mathematica, on a Pentium
Pro running under Linux.
In[1]:= normInfinity[m_] := Max[Apply[Plus, Abs[m], 1]]
In[2]:= norm1[m_] := normInfinity[Transpose[m]]
In[3]:= norm2[m_] := First[SingularValues[m][[2]]]
Let's try an example.
In[5]:= SeedRandom[1111]
In[6]:= mat = Table[Random[], {400}, {400}];
In[7]:= normInfinity[mat]
Out[7]= 216.233
In[8]:= norm1[mat]
Out[8]= 218.138
In[9]:= Timing[norm2[mat]]
Out[9]= {34.54 Second, 200.454}
Good so far, except that the 2-norm is time-consuming to compute. Below
is a faster way that tends to be quite accurate. We use the so-called
"power method" to find the largest eigenvalue of
mat.Transpose[Conjugate[mat]]. As this matrix is symmetric, its largest
eigenvalue is also its largest singular value. Moreover, this is just
the square of the largest singular value of mat.
In[10]:= estimateNorm2[mat_?MatrixQ]:= Module[
{symmat=mat.Transpose[Conjugate[mat]], len=Length[mat],
prec=Precision[mat], x, y},
y = Table[Random[Real, {1, 2}, prec+5], {len}];
For [j=1, j<=5, j++,
x = y;
y = symmat.x;
];
Sqrt[Apply[Plus, Map[Abs,y]] / Apply[Plus, Map[Abs,x]]]
]
In[11]:= Timing[estimateNorm2[mat]]
Out[11]= {13.11 Second, 200.454}
It can be checked that most of the time to do this is in taking the dot
product. Which makes sense, because that operation is O(n^3) for an
n-by-n matrix, whereas the rest of the code is clearly O(n^2).
Daniel Lichtblau
Wolfram Research
danl at wolfram.com