27
4
This challenge is to write fast code that can perform a computationally difficult infinite sum.
Input
An n
by n
matrix P
with integer entries that are smaller than 100
in absolute value. When testing I am happy to provide input to your code in any sensible format your code wants. The default will be one line per row of the matrix, space separated and provided on standard input.
P
will be positive definite which implies it will always be symmetric. Other than that you don't really need to know what positive definite means to answer the challenge. It does however mean that there will actually be an answer to the sum defined below.
You do however need to know what a matrix-vector product is.
Output
Your code should compute the infinite sum:
to within plus or minus 0.0001 of the correct answer. Here Z
is the set of integers and so Z^n
is all possible vectors with n
integer elements and e
is the famous mathematical constant that approximately equals 2.71828. Note that the value in the exponent is simply a number. See below for an explicit example.
How does this relate to the Riemann Theta function?
In the notation of this paper on approximating the Riemann Theta function we are trying to compute . Our problem is a special case for at least two reasons.
- We set the initial parameter called
z
in the linked paper to 0. - We create the matrix
P
in such a way that the minimum size of an eigenvalue is1
. (See below for how the matrix is created.)
Examples
P = [[ 5., 2., 0., 0.],
[ 2., 5., 2., -2.],
[ 0., 2., 5., 0.],
[ 0., -2., 0., 5.]]
Output: 1.07551411208
In more detail, let us see just one term in the sum for this P. Take for example just one term in the sum:
and x^T P x = 30
. Notice thate^(-30)
is about 10^(-14)
and so is unlikely to be important for getting the correct answer up to the given tolerance. Recall that the infinite sum will actually use every possible vector of length 4 where the elements are integers. I just picked one to give an explicit example.
P = [[ 5., 2., 2., 2.],
[ 2., 5., 4., 4.],
[ 2., 4., 5., 4.],
[ 2., 4., 4., 5.]]
Output = 1.91841190706
P = [[ 6., -3., 3., -3., 3.],
[-3., 6., -5., 5., -5.],
[ 3., -5., 6., -5., 5.],
[-3., 5., -5., 6., -5.],
[ 3., -5., 5., -5., 6.]]
Output = 2.87091065342
P = [[6., -1., -3., 1., 3., -1., -3., 1., 3.],
[-1., 6., -1., -5., 1., 5., -1., -5., 1.],
[-3., -1., 6., 1., -5., -1., 5., 1., -5.],
[1., -5., 1., 6., -1., -5., 1., 5., -1.],
[3., 1., -5., -1., 6., 1., -5., -1., 5.],
[-1., 5., -1., -5., 1., 6., -1., -5., 1.],
[-3., -1., 5., 1., -5., -1., 6., 1., -5.],
[1., -5., 1., 5., -1., -5., 1., 6., -1.],
[3., 1., -5., -1., 5., 1., -5., -1., 6.]]
Output: 8.1443647932
P = [[ 7., 2., 0., 0., 6., 2., 0., 0., 6.],
[ 2., 7., 0., 0., 2., 6., 0., 0., 2.],
[ 0., 0., 7., -2., 0., 0., 6., -2., 0.],
[ 0., 0., -2., 7., 0., 0., -2., 6., 0.],
[ 6., 2., 0., 0., 7., 2., 0., 0., 6.],
[ 2., 6., 0., 0., 2., 7., 0., 0., 2.],
[ 0., 0., 6., -2., 0., 0., 7., -2., 0.],
[ 0., 0., -2., 6., 0., 0., -2., 7., 0.],
[ 6., 2., 0., 0., 6., 2., 0., 0., 7.]]
Output = 3.80639191181
Score
I will test your code on randomly chosen matrices P of increasing size.
Your score is simply the largest n
for which I get a correct answer in less than 30 seconds when averaged over 5 runs with randomly chosen matrices P
of that size.
What about a tie?
If there is a tie, the winner will be the one whose code runs fastest averaged over 5 runs. In the event that those times are also equal, the winner is the first answer.
How will the random input be created?
- Let M be a random m by n matrix with m<=n and entries which are -1 or 1. In Python/numpy
M = np.random.choice([0,1], size = (m,n))*2-1
. In practice I will setm
to be aboutn/2
. - Let P be the identity matrix + M^T M. In Python/numpy
P =np.identity(n)+np.dot(M.T,M)
. We are now guaranteed thatP
is positive definite and the entries are in a suitable range.
Note that this means that all eigenvalues of P are at least 1, making the problem potentially easier than the general problem of approximating the Riemann Theta function.
Languages and libraries
You can use any language or library you like. However, for the purposes of scoring I will run your code on my machine so please provide clear instructions for how to run it on Ubuntu.
My Machine The timings will be run on my machine. This is a standard Ubuntu install on an 8GB AMD FX-8350 Eight-Core Processor. This also means I need to be able to run your code.
Leading answers
n = 47
in C++ by Ton Hospeln = 8
in Python by Maltysen
It may be worth mentioning that a positive definite matrix is symmetric by definition. – 2012rcampion – 2016-04-04T13:03:07.823
@2012rcampion Thanks. Added. – None – 2016-04-04T13:06:03.730
Ok, maybe this is a dumb question, but I've stared at this for ages and I cannot figure out how you got an
x
of[-1,0,2,1]
. Can you elaborate on this? (Hint: I am not a math guru) – wnnmaw – 2016-04-04T19:07:59.237@wnnmaw Sorry for being confusing. The sum has one term for every possible vector x of length 4 in this case. [-1,0,2,1] is just one I picked at random to show explicitly what the term would be in that case. – None – 2016-04-04T19:09:26.377
You're testing on randomly chosen matrices. A crucial thing to know: What distribution do they follow/how do you generate those? – flawr – 2016-04-04T22:27:20.673
Does the sum always converge? – soktinpk – 2016-04-05T00:22:04.707
@soktinpk Yes, because of the positive definite property. – None – 2016-04-05T06:41:13.900
@R.Kap M^T means "M transpose". https://chortle.ccsu.edu/VectorLessons/vmch13/vmch13_14.html
– None – 2016-04-05T08:30:45.597@R.Kap Yes. Vector transpose just swaps a vector between being a row vector and a column vector (or vice versa) – None – 2016-04-05T08:36:10.003
@R.Kap You can definitely use numpy. In fact I just removed the rule as it was a little confusing and it isn't needed. – None – 2016-04-05T08:45:32.103
@R.Kap The
e
denotes the exponential function. – flawr – 2016-04-05T17:55:26.080So the equation uses ALL possible vectors consisting of ALL the integers (whether positive or negative) with length
n
. Wow... – R. Kap – 2016-04-05T18:32:13.1201@Lembik The way you generate the SPD matrices implies that no singular value has ever an absolute value below 1. Can we use that knowledge? – flawr – 2016-04-05T18:43:47.450
@flawr Yes absolutely. – None – 2016-04-05T18:47:00.717
this one asks for matlab... gotta fire my instance up to set the arithmetic later – masterX244 – 2016-04-06T17:40:49.133
@masterX244 Please feel free to answer in matlab! I can't run matlab myself but I will try to run it in octave for scoring and even convert it to python if I can. That is not normally too hard these days. – None – 2016-04-06T17:42:12.977
So how do you actually know what the correct answer is? I mean, based on the equation given, from my perspective it seems as if their could be many different answers based on the bounds used. So how do you actually know what the correct answer (or range) should be? – R. Kap – 2016-04-09T08:44:50.687
@R.Kap Ton Hospel's code uses some provable bounds. Hopefully there will be more documentation soon. – None – 2016-04-09T09:40:24.987