Ponder This

May 2023 - Solution

<< April May June >>

May 2023 Challenge

May 2023 Solution:

The numerical results are
52390.274 for n=2^{20}
156032.269 for n=2^{30}
We have accepted other, less exact solutions.

The naive solution for computing x^TAx=\sum_{i,j=0}^{n-1}A_{ij}x_ix_j requires O(n^2) steps (ignoring k). Instead, We describe here a clever way to compute the form in O(n\cdot 2^k) steps, based on the solution communicated to us by Lorenzo Gianferrari Pini.

First of all, note that A_{ij} takes values from the length 2^k vector a, so we can rewrite the sum by agglomerating over the values of a:

\sum_{i,j=0}^{n-1}A_{ij}x_ix_j=\sum_{t=0}^{2^k-1}a_t\cdot\sum_{i,j:[Q_i==Q_j]=t}x_ix_j

Define a length 2^k vector u by u_{t}=\sum_{i,j:[Q_i==Q_j]=t}x_ix_j. Given u, the requested solution is a\cdot u^T .

From now on we can focus on computing u_{t} for a specific value of t. This is done using a technique similar to the inclusion-exclusion principle.

We illustrate this for k=5 and

In this case, \sum_{i,j:\left[Q_{i}==Q_{j}\right]=\left\{ 1,3,5\right\} }x_{i}x_{j}=x_{1}x_{2}+x_{2}x_{1} since Q_2 agree precisely on the indices 1,3,5 and disagree on 2,4. Since Q_3 Agrees with Q_1 on index 4 and with Q_2 on index 2, it does not participate in this sum.

We now present a different way of obtaining this sum (more efficient when used in general). For any subset A\subseteq\{1,\ldots,k\} Define an indicator function \left[Q_{i}=Q_{j}\right]^A taking the value 1 if Q_i and Q_j agree on the indices in A (indices not in A are not considered at all) and 0 otherwise.

Now note that the sum \sum_{i,j:\left[Q_{i}==Q_{j}\right]=\left\{ 1,3,5\right\} }x_{i}x_{j} can be written using indicators:

\sum_{i,j}x_{i}x_{j}\left[Q_{i}=Q_{j}\right]^{\left\{ 1\right\} }\left(1-\left[Q_{i}=Q_{j}\right]^{\left\{ 2\right\} }\right)\left[Q_{i}=Q_{j}\right]^{\left\{ 3\right\} }\left(1-\left[Q_{i}=Q_{j}\right]^{\left\{ 4\right\} }\right)\left[Q_{i}=Q_{j}\right]^{\left\{ 5\right\} }

Opening the product, we obtain a sum of products of the form

\sum_{i,j}x_{i}x_{j}\left(\left[Q_{i}=Q_{j}\right]^{\left\{ 1,3,5\right\} }-\left[Q_{i}=Q_{j}\right]^{\left\{ 1,2,3,5\right\} }- \left[Q_{i}=Q_{j}\right]^{\left\{ 1,3,4,5\right\} }+\left[Q_{i}=Q_{j}\right]^{\left\{ 1,2,3,4,5\right\} }\right)

Each summand is an indicator for a set A which is a superset of \{1,3,5\}. The sign of the summand is related to the oddity of the extra number of indices added to A (indices not from \{1,3,5\}, which come from a term of the form 1-\left[Q_{i}=Q_{j}\right]^{\left\{ r\right\} }). We can write this in short form as

\sum_{i,j}x_{i}x_{j}\sum_{\left\{ 1,3,5\right\} \subseteq A}\left(-1\right)^{\left|A\right|-\left|\left\{ 1,3,5\right\} \right|}\left[Q_{i}=Q_{j}\right]^{A}

Switching the summation order we obtain

\sum_{\left\{ 1,3,5\right\} \subseteq A}\left(-1\right)^{\left|A\right|-\left|\left\{ 1,3,5\right\} \right|}\sum_{i,j}\left[Q_{i}=Q_{j}\right]^{A}x_{i}x_{j}

We are left with the question of how to compute \sum_{i,j}\left[Q_{i}=Q_{j}\right]^{A}x_{i}x_{j}. This is relatively easy since we do not require Q_i, Q_j to be distinct in the indices not present in A. This means the relation i\sim_A j\iff\left[Q_{i}=Q_{j}\right]^{A}=1 is an equivalence relation over the set [n]=\{1,2,\ldots,n\}, and we can write:\sum_{i,j}\left[Q_{i}=Q_{j}\right]^{A}x_{i}x_{j}=\sum_{X\in\left(\left[n\right]/\sim_{A}\right)}\left(\sum_{i\in X}x_{i}\right)^{2}

In our example, we have
\sum_{i,j}\left[Q_{i}=Q_{j}\right]^{\left\{ 1,3,5\right\} }x_{i}x_{j}=\left(x_{1}+x_{2}+x_{3}\right)^{2} \sum_{i,j}\left[Q_{i}=Q_{j}\right]^{\left\{ 1,2,3,5\right\} } x_{i}x_{j}=x_{1}^{2}+\left(x_{2}+x_{3}\right)^{2} \sum_{i,j}\left[Q_{i}=Q_{j}\right]^{\left\{ 1,3,4,5\right\} } x_{i}x_{j}=\left(x_{1}+x_{3}\right)^{2}+x_{2}^{2} \sum_{i,j}\left[Q_{i}=Q_{j}\right]^{\left\{ 1,2,3,4,5\right\} } x_{i}x_{j}=x_{1}^{2}+x_{2}^{2}+x_{3}^{2}

And indeed
\left(x_{1}+x_{2}+x_{3}\right)^{2}-\left[x_{1}^{2}+\left(x_{2}+x_{3}\right)^{2}\right]- \left[\left(x_{1}+x_{3}\right)^{2}+x_{2}^{2}\right]+\left[x_{1}^{2}+x_{2}^{2}+x_{3}^{2}\right]=x_{1}x_{2}+x_{2}x_{1}

This works in general. For any t\in\left\{ 0,1,\ldots,2^{k}-1\right\} , denote by B_t the set of indices in the binary representation of t which are 1. We have

\sum_{i,j:\left[Q_{i}=Q_{j}\right]=t}x_{i}x_{j}=\sum_{B_{t}\subseteq A}\left(-1\right)^{\left(\left|A\right|- \left|B_{t}\right|\right)}\sum_{X\in\left(\left[n\right]/\sim_{A}\right)}\left(\sum_{i\in X}x_{i}\right)^{2}

Since the expression \sum_{X\in\left(\left[n\right]/\sim_{A}\right)}\left(\sum_{i\in X}x_{i}\right)^{2} is independent of B_t, it suffices to compute it once, for every subset A\subseteq\left\{ 0,1,\ldots,2^{k}-1\right\}. We create a vector v such that

v_{t}=\sum_{X\in\left(\left[n\right]/\sim_{A_{t}}\right)}\left(\sum_{i\in X}x_{i}\right)^{2}

The inclusion-exclusion part is handled by defining a matrix

M_{t,s}=\begin{cases}
\left(-1\right)^{\left|A_{s}\right|-\left|B_{t}\right|} & B_{t}\subseteq A_{s}\\
0 & B_{t}\not\subseteq A_{s}
\end{cases}

Recall that we defined
u_{t}=\sum_{i,j:[Q_i==Q_j]=t}x_ix_j
Such that a\cdot u^T is the required solution. We have now seen that
u=M\cdot v^T
Thus, the algorithm for solving the problem efficiently is

1. Compute the matrix M using the formula M_{t,s}=\begin{cases}
\left(-1\right)^{\left|A_{s}\right|-\left|B_{t}\right|} & B_{t}\subseteq A_{s}\\
0 & B_{t}\not\subseteq A_{s}
\end{cases}

2. Compute the vector v using the formula v_{t}=\sum_{X\in\left(\left[n\right]/\sim_{A_{t}}\right)}\left(\sum_{i\in X}x_{i}\right)^{2}

3. Compute a\cdot M\cdot v^T

Steps 1 and 3 involve operations polynomial in O(2^k); the main computation is done in step 2, where for each of the 2^k values of t, one has to go over the O(n) vector x and sum according to the equivalance classes.

An unoptimized Python code for this computation is

 
import functools 
import numpy as np 
k = 5 

@functools.lru_cache(maxsize=None) 
def Q(i): 
    return [int((2**k)* (np.sin((i+1)*(t+1))-np.floor(np.sin((i+1)*(t+1))))) for t in range(k)] 

def num_to_set(t): 
    binary = bin(t)[2:] 
    return {i for i, bit in enumerate(binary[::-1]) if bit == '1'} 

def M(t, s): 
    A = num_to_set(s) 
    B = num_to_set(t) 
    if not A.issuperset(B): 
        return 0 
    else: 
        diff = len(A.difference(B)) 
        if diff % 2 == 0: 
            return 1 
        else: 
            return -1 

def equiv_class_sum(t, x): 
    A = num_to_set(t) 
    results = {} 
    for i in range(n): 
        Qi = Q(i) 
        subQ = tuple([Qi[j] for j in A]) 
        if subQ not in results: 
            results[subQ] = 0 
        results[subQ] += x[i] 
    return sum(res**2 for res in results.values()) 

n = 2**30 
a = np.linspace(0, 1, 2**k).reshape((1,2**k)) 
M = np.array([[M(t,s) for s in range(2**k)] for t in range(2**k)]) 
x = np.linspace(-1,1,n) 
v = np.array([equiv_class_sum(t, x) for t in range(2**k)]).reshape((2**k,1)) 
print((a @ M @ v)[0][0])