9
$\begingroup$

Can you think of a way of speeding up this computation ? I don't.

n = 100000;
f[t1_, t2_] = 
  x^5*(x*(x*(x + I) + I) + I*(Exp[I*t2]*(Exp[I*t1]*(Exp[I*t2]*(-Exp[I*t2] + 1) + 1) - 1) - 
         1)) - I*x*(x*(x*(x - 1) - 1) + 1) + I*(Exp[I*t1]*(Exp[I*t1]*(Exp[I*t1]*(Exp[I*t2] - 1) - 1) + 1) - 1);
tlist = RandomReal[{0, 2 Pi}, {n, 2}];

data = Flatten[ParallelTable[ReIm[x /. NSolve[f @@ t == 0, {x}]], {t, tlist}], 1];

Why? To produce nice images!

ListPlot[data, PlotStyle -> {Opacity[0.025], PointSize[0.001], Black},
  AspectRatio -> Automatic, Axes -> False]

enter image description here

Reference: https://profconradi.github.io/category/matematica.html

$\endgroup$
9
  • $\begingroup$ How long does it take now? $\endgroup$ Commented Nov 27, 2023 at 19:11
  • 2
    $\begingroup$ You're basically finding the roots of a polynomial in x, so you can use Solve rather than NSolve. Doing so knocks about 15% off of the computer time on my machine; not a huge speed-up but not completely insignificant either. $\endgroup$ Commented Nov 27, 2023 at 19:15
  • 3
    $\begingroup$ This seems somewhat faster. nsols = ParallelTable[{ToRules[NRoots[f @@ t == 0, x]]}, {t, tlist}];. I get a bit under 3 seconds on my machine. $\endgroup$ Commented Nov 28, 2023 at 0:47
  • 2
    $\begingroup$ Strongly related: mathematica.stackexchange.com/q/73312/1871 $\endgroup$
    – xzczd
    Commented Dec 4, 2023 at 0:19
  • 1
    $\begingroup$ @PerAlexandersson The polynomial depends on parameters $(t_1,t_2)$. Each point corresponds to the roots of a polynomial for some random values of the parameters. $\endgroup$
    – anderstood
    Commented Dec 4, 2023 at 20:02

3 Answers 3

8
$\begingroup$

Since this is about polynomials in one complex variable, we can compute the eigenvalues of the companion matrix instead. Since Eigenvalues is a bit slow on small matrices, I decided to use a LibraryFunction. For ease of use, I employ the Eigen package for the eigensolve.

coeffs = Compile[{{t, _Real, 1}},
   With[{
     T1 = Exp[I Compile`GetElement[t, 1]],
     T2 = Exp[I Compile`GetElement[t, 2]]
     },
    {I (-1. + T1 (1. + T1 (-1. + T1 (-1. + T2)))), -I 1., I 1., 
     1. I, -1. I, 1. I (-1. + T2 (-1. + T1 (1. + (1. - T2) T2))), 
     1. I, 1. I, 1.}
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable}
   ];

Needs["CCompilerDriver`"];
Quiet[LibraryFunctionUnload[cSolvePolynomial[8]]];
ClearAll[cSolvePolynomial];
cSolvePolynomial[d_] := 
  cSolvePolynomial[d] = Module[{lib, code, name}, name = "cf";
    code = StringJoin["

#include \"WolframLibrary.h\"

#include <thread>
#include <future>
#include <Eigen/Eigenvalues>

static constexpr mint d = " <> IntegerString[d] <> ";

using Complex = std::complex<mreal>;

using Matrix_T = Eigen::Matrix<Complex,d,d>;
using Vector_T = Eigen::Matrix<Complex,d,1>;

    // Computes k-th job pointer for job_count equally sized jobs distributed on thread_count threads.
    template<typename Int, typename Int1, typename Int2>
    inline Int JobPointer( const Int job_count, const Int1 thread_count_, const Int2 thread_ )
    {
        const Int thread_count = static_cast<Int>(thread_count_);
        const Int thread       = static_cast<Int>(thread_);
        const Int chunk_size   = job_count / thread_count;
        const Int remainder    = job_count % thread_count;
        
        return chunk_size * thread + (remainder * thread) / thread_count;
    }


    // Executes the function `fun` of the form `[]( const Int thread ) -> void {...}` parallelized over `thread_count` threads.
    template<typename F, typename Int>
    inline void ParallelDo( F && fun, const Int thread_count )
    {
        if( thread_count <= static_cast<Int>(1) )
        {
            std::invoke( fun, static_cast<Int>(0) );
        }
        else
        {
            std::vector<std::future<void>> futures (thread_count);
            
            for( Int thread = 0; thread < thread_count; ++thread )
            {
                futures[thread] = std::async( std::forward<F>(fun), thread );
            }
            
            for( auto & future : futures )
            {
                future.get();
            }
        }
    }

EXTERN_C DLLEXPORT int " <> name <> 
       "(WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res)
{
    MTensor coeffs_      = MArgument_getMTensor(Args[0]);
    mint    thread_count = MArgument_getInteger(Args[1]);

    const mint n = libData->MTensor_getDimensions(coeffs_)[0];

    const Complex * const coeffs = reinterpret_cast<Complex*>(libData->MTensor_getComplexData(coeffs_));

    MTensor roots_;

    // Instead of a complex array of size n x d we allocate a real array of size (n * d) x 2.
    // This way, we can get rid of ReIm and Flatten.

    const mint dims [2] { n * d, 2 };

    (void)libData->MTensor_new( MType_Real, 2, &dims[0], &roots_ );

    Complex * roots = reinterpret_cast<Complex*>( libData->MTensor_getRealData(roots_));

    // This loops over all coefficient vectors and computes the eigenvalues of the companion matrix.

    ParallelDo(
        [roots,coeffs,n,thread_count]( const mint thread )
        {
            Eigen::ComplexEigenSolver<Matrix_T> S;
    
            Matrix_T C;
    
            C.setZero();
        
            for( int i = 0; i < d-1; ++i )
            {
                C(i+1,i) = Complex { 1., 0. };
            }
    
            const mint k_begin = JobPointer( n, thread_count, thread     );
            const mint k_end   = JobPointer( n, thread_count, thread + 1 );
    
            for( mint k = k_begin; k < k_end; ++k )
            {
                const Complex * const c = &coeffs[(d+1) * k];
                      Complex * const r = &roots [ d    * k];
        
                const Complex factor = - 1. / c[d];
        
                for( int i = 0; i < d; ++i )
                {
                    C(i,d-1) = c[i] * factor;
                }
        
                S.compute(C);
        
                Vector_T eigs = S.eigenvalues().col(0);
        
                for( int i = 0; i < d; ++i )
                {
                    r[i] = eigs(i);
                }
            }
        },
        thread_count
    );

 
    MArgument_setMTensor(Res, roots_);

    return LIBRARY_NO_ERROR;
}"];
    lib = 
     CreateLibrary[code, name, "Language" -> "C++", 
      "CompileOptions" -> {"-O3", "-std=c++20", "-pthread"}, 
      "IncludeDirectories" -> {"/opt/homebrew/include/eigen3"}, 
      "ShellOutputFunction" -> (If[# =!= "", Print[#]] &)];
    LibraryFunctionLoad[lib, 
     name, {{Complex, 2, "Constant"}, Integer}, {Real, 2}]
    ];

Here an usage example :

degree = 8;
threadCount = 8;
cf = cSolvePolynomial[degree];
data2 = cf[coeffs[tlist], threadCount]; // AbsoluteTiming // First

0.238111

On my machine, OP's code takes 52.2739 seconds instead. So this is a 220-fold speedup.

Correctness can be checked with

Max[Abs[(Sort /@ Partition[data, 8]) - (Sort /@ Partition[data2, 8])]]

1.05471*10^-14

This allows me to run this with ten times the number of points and to get this pictures:

enter image description here

And here with 8000000 random points:

enter image description here

$\endgroup$
6
$\begingroup$
data1 = Flatten[Table[ReIm[x /. NSolve[f @@ t == 0, {x}]], {t, tlist}], 1]; // AbsoluteTiming // First
(*    51.9531    *)

Optimizing @DanielLichtblau's comment a bit, and using $s_k=e^{i t_k}$:

g[s1_, s2_] = Series[f[t1, t2] /. {E^(I t1) -> s1, E^(I t2) -> s2},
                {x, 0, 8}] // FullSimplify // Normal
(*    I (-1 + s1 + s1^2 (-1 + s1 (-1 + s2))) - I x + I x^2 + I x^3 - I x^4 +
      I (-1 + s2 (-1 + s1 (1 + s2 - s2^2))) x^5 + I x^6 + I x^7 + x^8    *)

data2 = ReIm[List @@ Join @@ Map[NRoots[g @@ #, x] &, Exp[I*tlist]][[All, All, 2]]]; // AbsoluteTiming // First
(*    3.71658    *)

is faster by a factor of 14, and numerically indistinguishable:

data2 - data1 // Norm
(*    5.06198*10^-14    *)

Replacing Map by ParallelMap gives a small speedup, but not much.

$\endgroup$
5
$\begingroup$

Another possibility is the Ehrlich-Aberth algorithm which is easy to implement:

cAberth = Block[{z, t, x},
  With[{
    \[Epsilon] = 16. $MachineEpsilon,
    d = 8,
    rCode = N[
       Divide[
         f[Compile`GetElement[t, 1], Compile`GetElement[t, 2]], 
         D[f[Compile`GetElement[t, 1], Compile`GetElement[t, 2]], x]
         ] /. x -> z]
    },
   Compile[{{t, _Real, 1}},
    Block[{z, w, sum, maxw, r, zlist, wlist},
     zlist = RandomComplex[{-1 - I, 1 + I}, d];
     z = 0. + 0. I;
     r = 0. + 0. I;
     maxw = 1.;
     
     While[ maxw > \[Epsilon],
      
      maxw = 0.;
      
      Do[
       z = Compile`GetElement[zlist, i];
       r = rCode;
       sum = 0. + 0. I;
       Do[sum += 1./(z - Compile`GetElement[zlist, j]), {j, 1, i - 1}];
       Do[sum += 1./(z - Compile`GetElement[zlist, j]), {j, i + 1, d}];
       w = r /(1. - r sum);
       maxw = Max[maxw, Abs[w]];
       zlist[[i]] -= w;
       
       , {i, 1, d}];
      ];
     zlist
     ],
    CompilationTarget -> "C",
    RuntimeAttributes -> {Listable},
    Parallelization -> True,
    RuntimeOptions -> "Speed"
    ]
   ]
  ];

Use it like this on you problem:

data3 = Flatten[ReIm[cAberth[tlist]], 1];

This is a quick and dirty implementation of the "Gauss-Seidel"-like variant as mentioned in the text. I just generate random starting values without looking at the actual coefficients; so this should definitely be improved.

It takes about 25% longer than SolvePolynomial[degree][coeffs[tlist], threadCount] from my other answer. But I think that is due to the parallelization in CompiledFunctions; that is sometimes a bit poor. I expect that a corresponding LibraryLink implementation with parallelization within the function would perform better.

$\endgroup$
1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.