       Re: Re: Product command with matrices

• To: mathgroup at smc.vnet.net
• Subject: [mg87906] Re: [mg87889] Re: Product command with matrices
• From: danl at wolfram.com
• Date: Sat, 19 Apr 2008 23:50:41 -0400 (EDT)
• References: <fu9ft7\$ce2\$1@smc.vnet.net> <fu9vr5\$iim\$1@smc.vnet.net>

```> On Apr 18, 6:14 am, Jens-Peer Kuska <ku... at informatik.uni-leipzig.de>
> wrote:
>> Hi,
>>
>> Times[] is not Dot[] and Product[] and no equivalent to use
>> Dot[] instedad of Times[] and there is no (easy) way to
>> tell Mathematica that f is a matrix and what may be A
>> a scalar ? a vector or a matrix too ? Even I can't find it out
>> and Mathematica can't know it.
>>
>> You mean
>>
>> DotProduct[mtx_, {i_, i1_, in_}] := Dot @@ Table[mtx, {i, i1, in}]
>>
>> Phi := Id;
>> Phi[n_] := DotProduct[Id + f[i] ** A, {i, 0, n - 1}]
>>
>> and
>>
>> Phi gives
>>
>> (Id + f ** A).(Id + f ** A).(Id + f ** A)
>>
>> Regards
>>    Jens
>>
>> J Davis wrote:
>> > I want to define a matrix valued function such as the following (in
>> > LaTeX lingo):
>>
>> > \$\$
>> > X(0)=Id,
>> > X(n)=\prod_{i=0}^{n-1} (Id + f[i] A)
>> > \$\$
>>
>> > where A and f have already been defined and Id is the identity matrix
>> > of appropriate size.
>>
>> > I tried the following:
>>
>> > Id=IdentityMatrix;
>> > Phi:=Id;
>> > Phi[n_]:= Product[Id + f[i] A,{i,0,n-1}]
>>
>> > However, Phi and (Id+fA).(Id+fA).(Id+fA) do not agree.
>>
>> > Any help around this would be appreciated.
>
> Thanks very much to W Craig Carter (and others) who replied to me via
> email. I thought I would share the following very efficient solution:
>
> Phi[n_]:=Apply[Dot,Table[Id + f[k]A,{k,0,n-1}]]

There is a caveat to this approach. If your purpose is to do
matrix-times-vector products (this is quite common, say in iterative
linear solvers), then you'd want to avoid the explicit Dot expansion. That
is, instead of the matrix-times-matrix you would only do
matrix-times-vector. This will cut the run time by a factor approximating
the dimension. Here is an example, using integers just to show equivalence
(for numeric purposes, say, one would more likely have matrices of inexact
numbers).

f[j_] := 3*j^2 - j + 5
matfunc[m_, n_] := Table[IdentityMatrix[dim] + f[j]*m, {j, n}]

In:= dim = 25;
n = 200;
bnd = 5;
mat = RandomInteger[{-bnd, bnd}, {dim, dim}];
vec = RandomInteger[{-bnd, bnd}, {dim}];

Now notice that creating the table is quick, and evaluating the
matrix-times-vector is nearly dim times slower than expanding the matrix
product explicitly followed by application to the vector (and they give
the same result).

In:= Timing[mtable = matfunc[mat, n];]
Timing[prod = Apply[Dot, mtable];]
Timing[newvec = prod.vec;]
Timing[newvec2 = Fold[Dot[#2, #1] &, vec, Reverse[mtable]];]
Max[Abs[newvec - newvec2]]

Out= {0.02, Null}

Out= {3.785, Null}

Out= {0., Null}

Out= {0.171, Null}

Out= 0

Daniel Lichtblau
Wolfram Research

```

• Prev by Date: Re: axis alignment of 3D plots with ListContourPlot3D
• Next by Date: Re: Re: A kernel, multiple notebooks, and Global?
• Previous by thread: Re: Product command with matrices
• Next by thread: Scaling Plot inside Row