Benchmark Computations
A.I. McLeod

Cargo Loading Problem

The problem is to find the optimal loading of a carrier with n different items. Each item has value (Subscript[v, i]) and weight (Subscript[w, i]). The total weight capacity of the carrier is k. In addition, we suppose that (Subscript[q, i])  ≥ 1 items of type i are available. So as a mathematical programming problem we can state this as follows: find the optimal values for the non-negative integer decision variables, (Subscript[x, 1]) ,   ...  ,  (Subscript[x, n])  such that

max (∑_(i = 1)^n) (Subscript[v, i]) (Subscript[x, i])

subject to the constraints

(∑_(i = 1)^n) (Subscript[w, i]) (Subscript[x, i])  ≤ k

and

0 ≤ x ≤  (Subscript[q, i . ])

The following Mathematica algorithm uses recursion to solve this problem. It is worth noting that the computation is greatly speeded up by saving previously computed function values. Many similar recursive procedures can be made more efficient using this technique.

\(Cargo[k_Integer?Positive, \ v_List?VectorQ, \ w_List?VectorQ, \ 
        q_List?VectorQ]\  := 
\ 
      Module[{}, 
        
 (*\ check\ that\ all\ vectors\ are\ the\ same\ length\ *) 
\ \ 
        If[Map[Length, \ {v, w, q}] === RotateLeft[Map[Length, \ {v, w, q}]], 
          
\ \ \ \ \ \ Null, 
          Print["Error: vector arguments lengths unequal"]; 
          Return[Null]\ ]; 
        
 (*\ check\ that\ all\ coefficients\ are\ positive\ *) 
\ \ 
        If[Apply[And, \ Flatten[Positive[{v, w, q}]]], Null, 
\ \ \ \ \ 
          Print["Error: coefficient negative or zero"]; Return[Null]\ ]; 
        
 (*\ for\ speed, \ v2, \ w2, \ 
          q2\ are\ passed\ as\ static\ variables\ to\ f\ *) 
\ \ v2\  = \ v; 
        
\ \ w2\  = \ w; 
\ \ q2\  = \ q; \ 
\ \ f[k, \ Length[v]]

        \ \ \ \ \ \ \ ]; 
f[k_, \ 0]\  = \ {0}; 

    f[k_, \ i_]\  := \ 
\ 
      Module[\ {z, \ t, \ opt}, 
\ \ \ \ \ \ \ 
        opt\  = \ Array[0&, \ {i + 1}]; 
\ \ \ \ \ \ \ 
        Do[opt = 
            If[
\ \ \ \ \ \ \ \ \ \ \ \ \ 
              \((z\  = \ 
                    x\ v2[\([i]\)]\  + \ 
                      First[t = f[k - w2[\([i]\)]\ x, \ i - 1]])\)\  > \ 
                First[opt], \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ 
              Join[{z}, Rest[t], {x}], \ opt], 

          \ \ \ \ \ \ \ \ \ \ \ \ \ \ {x, 0, 
            Min[q2[\([i]\)], \ Floor[k/w2[\([i]\)]]]}
\ \ \ \ \ \ \ \ \ ]; 
        
 (*\ note : \ 
            this\ assignment\ reduces\ cpu\ time\ by\ a\ factor\ of\ 6*) 

        \ \ \ \ \ \ \ Return[f[k, i] = opt]
\ \ \ \ \ \ \ \ \ ]; \)

\(v\  = \ {72, \ 60, \ 40, \ 27, \ 20, \ 50, \ 85, \ 96}; \ 

    w\  = \ {20, \ 18, \ 14, \ 12, \ 10, \ 16, \ 22, \ 24}; \ 

    q\  = \ {99, \ 99, \ 99, \ 99, \ 99, \ 99, \ 99, \ 99}; \)

Timing[Cargo[181, v, w, q]]

Kernel Density Estimation and Plotting

The general form of the kernel estimator may be written,

\(TraditionalForm
    \`\(\(\(f\&^\)(x)\  = \ 
      1\/\(n\ h\)\ \(∑\+\(i = 1\)\%n K(\(x - X\_i\)\/h)\), \)\ \)\)

where

\(TraditionalForm
    \`∫\_\(-∞\)\%∞\( K(t)\) 
          dt\  = \ 1, \ \ \ \ \ K(t)\  ≥ \ 0, 
    \ \ ∀ \ \(t . \)\)

The parameter h is the bandwidth. Sometimes it is also called the window width or smoothing parameter. A popular kernel is the standard Gaussian

\(TraditionalForm
    \`K(t)\  = \ \(1\/\@\(2\ π\)\) e\^\(\(-\(1\/2\)\) t\^2\)\)

When  K(t) is chosen to be a smooth bell-shaped curve such as the standard Gaussian or the Tukey bisquare, then smooth estimators of the probability density function are obtained. It has been shown that the exact choice of the shape K(t) is not very important so long as it is smooth and bell-shaped (Silverman, 1986). The choice of the smoothing parameter, h, is very important and there is no firm rule. One chooses h keeping in mind the trade-off between too much variability and roughness on the one hand and like of fidelity and increased bias on the other hand. This is the old story: bias vs. variance. The best choice of h requires some experimentation by the data analyst. Silverman (1986, p.47) has shown that a good initial choice or rule-of-thumb for use with the Gaussian kernel is given by (Subscript[h, opt]) =  0.9 A (n^(-1/5)), where   A  =  min( standard deviation, interquartile range/1.34 ).
Just as the naive estimator could be viewed as the sum of rectangular boxes, the kernel estimator with a standard Gaussian kernel can be viewed as the sum of Gaussian bumps. 
My Mathematica implementation of a density smoother uses the Gaussian kernel and default smoothing given by Silverman's rule-of-thumb.  We try it out on the rivers and returns data. Unfortuantely it takes a lot of time!!

Needs["Statistics`DescriptiveStatistics`"]

\(PDFPlot[data_, \ spar_: 1]\  := \ 
\ \ 
        Module[{x\  = \ N[data]}, \ 
          hopt\  = \ 
\ \ \ \ \ 
            \((0.9*Min[
                    \((Quantile[x, \ 0.75]\  - \ Quantile[x, \ 0.25])\)/1.34, 
                    \ 
\ \ \ \ \ \ \ \ \ Sqrt[Variance[x]]]*spar)\)/
              Length[x]^0.2; \ 
\ \ \ \ 
          dens[z_, \ h_, \ x_]\  := \ 
\ \ \ \ \ 
            Plus\ @@\ 
              \(\((N[1/Sqrt[2*Pi]]/E^\((0.5*#1^2)\)\ &\ )\)\ /@\ 

                \ \ \ \ \ \ \ 
                Select[Table[\((z\  - \ x[\([i]\)])\)/h, 
                    \ {i, \ 1, \ Length[x]}], \ #1\  < \ 2.5\ &\ ]\); \ 

          \ \ \ \ Plot[dens[z, \ hopt, \ x], \ {z, \ Min[x], \ Max[x]}, \ 
            PlotRange\  -> \ All]]\)

Timing[PDFPlot[Table[Random[] + Random[], {300}]]]

Factorization

Factorization of a long integer tests integer arithmetic operations.

Timing[FactorInteger[5^100 - 1]]

Some Typical Benchmarks

({{Test, Machine, Timing }, {Cargo, Thinkpad,  Pentium - 120, 4.45}, {PDFPlot,  Thinkpad,  Pentium - 120, 78.1}, {FactorInteger, Thinkpad,  Pentium - 120, 16.2}, {Cargo, Risk 600 M40, (lexis), 4.41}, {PDFPlot, Risk 600 M40, (lexis), 98.91}, {FactorInteger, Risk 600 M40 (lexis), 15.33}})