0
$\begingroup$

I will do my best to keep this question as short as possible, I am afraid though that I wont be able to get my issue across if I don't elaborate on the background and provide my code as an example.

Issue: I enjoy using Mathematica for my Numerics homework. I have a hard time using Matlab because I always need to translate things (i.e. algorithms and closed formulas) first into Matrices and Vectors before I can turn to the actual problem.

Mathematica instead lets me dive right into the implementation, the price I pay for it is the performance. I have no doubt that this is my fault hence I am here to ask for advice.


My issue as described above has shown dramatically when I was implementing the Discrete Fourier Transformation (of a real valued function) for a class assignment (thus not allowed to use built in functions). I was given these two items:

Def: Let $f_r$, $0 \leq r \leq N-1$ be the values of a periodic function with length $b-a$. $$ \hat{f_r}:= \frac{1}{N} \sum_{m=0}^{N-1} f_m e^{- \frac{ 2 \pi i}{N} m r}, \ 0 \leq r \leq N-1 $$


Lemma: Let $f:[a,b] \to \mathbb{R}$ with data $f_r, 0 \leq r \leq N-1$, then the trigonometric Fourier Transformation is given by $$T(t) = \hat{A_0} + \sum_{m=1}^{\left\lfloor{\frac{N-1}{2}} \right\rfloor} \left( \hat{A_m} \cos \left( \frac{2 \pi}{b-a}mt \right) + \hat{B_m} \sin \left( \frac{2 \pi }{b-a}mt \right) \right) \\ + \begin{cases} \hat{A_{N/2}} \cos \left( \frac{ \pi N}{b-a}t\right) \text{, if $N$ is even} \\ 0 \text{, if $N$ is odd}\end{cases} $$ Where $\hat{A_0}=\hat{f_0}, \ \hat{A_{N/2}}=\hat{f_{N/2}}$ and
$\hat{A_m} = 2 \text{Re}( \hat{f_m}), \ \hat{B_m} = -2 \text{Im}(\hat{f_m})$

Since these are beautifully closed formulas I would indeed implement them using Mathematica, naively following what the definition/lemma told me to do.

The issue I encounter is when $N$ becomes large, I have already massive trouble when $N=2^{12}=4096$

My implementation:

fhatentry[r_, rvector_] := 
 fhatentry[r, 
   rvector] = (1./Length[rvector]) Sum[ 
    rvector[[m + 1]] Exp[-(2.*Pi*I*m*r)/Length[rvector]], {m, 0., 
     Length[rvector] - 1}]

As given by the Definition, using rvector for $f_r$. Note that I have already tried improving my poor performance speed by forcing the function to recall its values (might be not a good idea) and using numerical values, i.e. instead of the symbol $2$, using it's numerical equivalent $2.$

As in the Lemma I have implemented

ahat[k_, rvector_] := 2 * Re[ fhatentry[k, rvector]]
bhat[k_, rvector_] := -2* Im[fhatentry[k, rvector]]

and

trigpol[t_, rvector_] := 
 fhatentry[0, rvector] + 
  Sum[ ahat[m, rvector] Cos[ (2. Pi m t)/(b - a)] + 
    bhat[m, rvector] Sin[ (2. Pi m t)/(b - a)], {m, 1., Floor[ (Length[rvector] - 1)/2]}] + 
  If[ EvenQ[Length[rvector]], 
   fhatentry[ Length[rvector]/2, 
     rvector] Cos[ (Pi Length[rvector] t)/(b - a)], 0]

I have also given that $a=0,b=1$ and stored them in Mathematica accordingly using numerical values.


Question: I am given a vector, say data, with $4096$ entries which translate to the values of $f_r$ for some periodic function $f:[0,1] \to \mathbb{R}$ and I am interested in the plot in the region [0,1].

Plot[trigpol[t,data], {t,0,1}]

Does terminate and does its job, but it takes like 5 minutes to do so.

I understand if these kind of questions are tedious but I am curious on how to:

  • Tune the performance of my code above to make it work faster, while still keeping it relatively legible and easy to understand what's going on.

Update

In the comments I was suggested to use a vectorized coding style to implement my functions, I have done that as follows

trigpol2[t_, rvector_] := With[{len = Length[rvector]}, 
  fhatentry[0, rvector] + 
   Total[ 
    ahat[#, rvector] & /@ Range[1, Floor[(len - 1)/2]] *
     Cos[(2 Pi Range[1, Floor[(len - 1)/2]] t )/(b - a)]] + 
   Total[ bhat[#, rvector] & /@ Range[1, Floor[(len - 1)/2]]* 
     Sin[( 2 Pi Range[1, Floor[(len - 1)/2]] t)/(b - a)]] 
   + If[ EvenQ[len], 
    fhatentry[ len/2, rvector] Cos[ (len t)/(b - a)], 0]
  ]

When plotting this function (as in the Question above) my Timing went from 198 to 251 (the output remains correct though)

$\endgroup$
  • $\begingroup$ try Plot[Evaluate[trigpol[t, data]], {t, 0, 1}] or Plot[trigpol[t, data], {t, 0, 1}, Evaluated -> True] $\endgroup$ – kglr Apr 25 '17 at 2:05
  • $\begingroup$ Vectorized operations on numerical data tends to be much faster than Summing element-wise, so we can try to implement this without losing too much readability. I would replace your fhatentry with fhatentry[r_, rvector_] := With[{len = Length[rvector]}, Total[rvector*Exp[-2 Pi I Range[0, len - 1]*r/len]]/len]. There is almost no point in memoizing this function because each value is only used twice: once for ahat and once for bhat. I'll let you think about how to "vectorize" the Sum in trigpol :) $\endgroup$ – Marius Ladegård Meyer Apr 25 '17 at 7:05
  • $\begingroup$ Dear @MariusLadegårdMeyer - Thanks a lot for your advice, I tried my best to follow it and I implemented the vectorized version for trigpol, I will post it in the OP accordingly. I am sure the fault is entirely on my side, but it didn't help reduce the computation time, in fact it increased it. $\endgroup$ – Spaced Apr 25 '17 at 13:38
  • $\begingroup$ Also, I think you will have more control over the time the plotting takes if you use ListPlot instead, where you choose the t points yourself. E.g. ListPlot[Table[{t,trigpol[t,data]},{t,0,1,1/99}]] will give you 100 plot points evenly spaced. Adjust accordingly. $\endgroup$ – Marius Ladegård Meyer Apr 25 '17 at 13:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.