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[1] 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[0] := Id; >> Phi[n_] := DotProduct[Id + f[i] ** A, {i, 0, n - 1}] >> >> and >> >> Phi[3] gives >> >> (Id + f[0] ** A).(Id + f[1] ** A).(Id + f[2] ** 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[2]; >> > Phi[0]:=Id; >> > Phi[n_]:= Product[Id + f[i] A,{i,0,n-1}] >> >> > However, Phi[3] and (Id+f[2]A).(Id+f[1]A).(Id+f[0]A) 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[82]:= 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[87]:= 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[87]= {0.02, Null} Out[88]= {3.785, Null} Out[89]= {0., Null} Out[90]= {0.171, Null} Out[91]= 0 Daniel Lichtblau Wolfram Research