1
$\begingroup$

edit2: The tweak suggested by Henrik Schumacher in Carl Woll's answer is so far a clear winner ($\times800$ over the benchmark); it's almost disturbing to contemplate how 'slow'-relatively speaking-something like -(list - center) is...

edit1:

So far this is the best alternative (without compiling anything yet)

negativeSemivariance5 = Function[{list, center, length}, 
  Total[Ramp[-(list - center)]^2]/(length - 1)
 ]

Evaluating the testing code snippet (simply replace negativeSemivariance4 with negativeSemivariance5) returns

{190., True}

which indicates a $\times190$ improvement on the specified test case.


original:

Is there a better way to write a semivariance function?

This is my benchmark definition

negativeSemivariance0[list_, center_, length_] := Module[{sum = 0.},
   Scan[If[# < center, sum += (# - center)^2] &, list];
   sum/length
  ]

and this is the best version I can think of

negativeSemivariance4[list_, center_, length_] := 
 With[{first = First[list]},
  Fold[
    With[{d = center - #2},
      #1 + Ramp[d]^2
      ] &,
    Ramp[center - first]^2,
    Rest[list]
    ]/length
  ]

I use the following code for testing:

Through[{Divide @@ Part[#, All, 1] &, Equal @@ Part[#, All, -1] &}[BlockRandom[
 With[{n = 5000},
  With[{rand = RandomReal[{-1, 1}, n]},
   With[{mean = Mean[rand]},
    {negativeSemivariance0[rand, mean, n] // RepeatedTiming, 
     negativeSemivariance4[rand, mean, n] // RepeatedTiming}
    ]
   ]
  ], RandomSeeding -> 132456987]]]

Evaluating it, returns

{20., True}

which can be interpreted as saying that the second version is $\times 20$ faster than the first and both functions return the same result.

Ps.This function will be frequently called on many different lists of same length. The size of the list used for testing-n=5000-is an average problem instance size.

$\endgroup$
2
$\begingroup$

I think using Dot should provide a speed boost:

svariance = Function[{list, center, length}, 
    #.#&[Ramp[-(list - center)]]/(length - 1)
];

Comparison:

data = RandomReal[{-1,1}, 10^6];

svariance[data, .1, 10^6-1] //RepeatedTiming
negativeSemivariance5[data, .1, 10^6-1] //RepeatedTiming

{0.0026, 0.222223}

{0.0050, 0.222223}

$\endgroup$
  • 1
    $\begingroup$ Function[{list, center, length}, #.# &[Ramp[Subtract[center, list]]]/(length - 1)] seems to be even a tick faster. $\endgroup$ – Henrik Schumacher Jun 3 '18 at 23:29
1
$\begingroup$
Variance[Select[data, # < Mean[data] &]]
$\endgroup$
  • 1
    $\begingroup$ Variance[] uses Mean[Select[]] which is not the same as Mean[data]; the denominator used for calculating variance is Length[Select[]]-1 instead of Length[data]-1 $\endgroup$ – user42582 Feb 12 '18 at 17:29
1
$\begingroup$

Replacing Ramp[-(list-center)] with Ramp[Subtract[center,list]] in @Carl's version seems to give a slight improvement in timings. Also, if the third argument length can be taken as the length of the input list, you can take the second Moment of the clipped list to get timings better than negativeSemivariance5.

semivar = Function[{list, center, length}, #.#&@Ramp[Subtract[center, list]]/(length-1)];
moment1 = Function[{list, center}, 
   MomentEvaluate[Moment[2], Ramp[Subtract[center, list]]]]; 
 moment2 = Function[{list, center}, Moment[Ramp[Subtract[center, list] ], 2]]; 

Timings:

SeedRandom[1]
data = RandomReal[{-1, 1}, 10^6]; 
{RepeatedTiming @ semivar[data, .1, 10^6 ] , 
 RepeatedTiming @ svariance[data, .1, 10^6 ] ,
 RepeatedTiming @ moment1[data, .1 ] ,  
 RepeatedTiming @ moment2[data, .1 ] ,   
 RepeatedTiming @ negativeSemivariance5[data, .1, 10^6 ]} 

{{0.0017, 0.222018}, {0.0021, 0.222018}, {0.055, 0.222017}, {0.0061, 0.222017},
{0.0066, 0.222018}}

$\endgroup$
  • $\begingroup$ yes indeed length is equal to the length of list; the tweak in Ramp is actually the one with the best performance so far; I was totally oblivious about the performance gains obtainable with Moment (needs a slight correction by a factor n/(n - 1); it performs slightly worse than negativeSemivariance5 but it definitely is an improvement over the benchmark) $\endgroup$ – user42582 Jun 4 '18 at 20:26

Your Answer

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

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