Search code examples
performancejuliaquadratic-programming

How can I efficiently calculate a quadratic form in Julia?

I would like to calculate a quadratic form: x' Q y in Julia.
What would be the most efficient way to calculate this for the cases:

  1. No assumption.
  2. Q is symmetric.
  3. x and y are the same (x = y).
  4. Both Q is symmetric and x = y.

I know Julia has dot(). But I wonder if it is faster than BLAS call.

Solution

  • Julia's LinearAlgebra stdlib has native implementations of 3-argument dot, and also a version that is specialized for symmetric/hermitian matrices. You can view the source here and here.

    You can confirm that they do not allocate using BenchmarkTools.@btime or BenchmarkTools.@ballocated (remember to interpolate variables using $). The symmetry of the matrix is exploited, but looking at the source, I don't see how x == y could enable any serious speedup, except perhaps saving a few array lookups.

    Edit: To compare the execution speed of the BLAS version and the native one, you can do

    1.7.0> using BenchmarkTools, LinearAlgebra

1.7.0> X = rand(100,100); y=rand(100);

1.7.0> @btime $y' * $X * $y
  42.800 μs (1 allocation: 896 bytes)
1213.5489200642382

1.7.0> @btime dot($y, $X, $y)
  1.540 μs (0 allocations: 0 bytes)
1213.548920064238

    This is a big win for the native version. For bigger matrices, the picture changes, though:

    1.7.0> X = rand(10000,10000); y=rand(10000);

1.7.0> @btime $y' * $X * $y
  33.790 ms (2 allocations: 78.17 KiB)
1.2507105095988091e7

1.7.0> @btime dot($y, $X, $y)
  44.849 ms (0 allocations: 0 bytes)
1.2507105095988117e7

    Possibly because BLAS uses threads, while dot is not multithreaded. There are also some floating point differences.