# How to model and solve such a 0-1 programming problem

My problem is described in this picture(It's like a Pyramid structure): The objective function is below: $$\min\sum_{k=1}^\ell\sum_{i=0}^{2^k-1}\sum_{j=0}^{2^k-1}\left(A_{i,j}^k-A_{i,j+1}^k\cdot\frac{D_{i,j}^k}{D_{i,j+1}^k}\right)^2$$

In this problem, $$D$$ is known, $$A$$ is the object that I want to get. It is a layered structure, each block in the upper layer is divided into four sub-blocks in the layer below. And the value of the upper layer node is equal to the sum of the four child nodes of the lower layer. In above example, I used only 2 layers.

What I want to do is simulate the distribution of $$D$$ with $$A$$, so in the objective function is the ratio of two adjacent squares in each row in $$A$$ compared to the value in $$D$$. I do this comparison on each layer and sum them. Then it is all of my objective function. But in the finest layer, the value in $$A$$ is $$0$$ or $$1$$ like the above picture, so this is also a 0-1 programming problem. I have tried to solve it using integer quadratic programming in python library CVXPY but adding the constraints of A<=1. However, it seems the speed is slow.

So I want to model this problem as a standard 0-1 programming problem, and then use some solutions of 0-1 programming to solve it, but I don't know how to model this problem as a standard 0-1 programming problem, and what methods I can use to solve it?

In other way, if I use quadratic programming instead of Integer quadratic programming (from $$A$$ is $$0/1$$ to $$A$$ is a value between $$0$$ and $$1$$). In this way whether I can solve it using the method of derivation, because this is a convex optimization problem, which can guarantee the global optimal solution. There are two unknown variables in each item, that is, the two items with $$A$$ in the formula. Partial derivatives are obtained for them, and the restriction of A<=1 is added, then solve using (such as) gradient descent method. Is this mathematically feasible, because I don't know much about optimization, and if it is possible, how should I do it? If not possible, what other methods can I use?

Say the finest layer, $$k$$ is 1, then declare $$n \times m$$ variables $$a_{n,m}^k \in\ \{0,1\}$$ as binary for the finest layer. For other layers, $$k=2...K$$
$$A_{i,j}^{k+1} = a_{i,j}^k+a_{i+1,j+1}^k+a_{i+1,j}^k+a_{i,j+1}^k$$.
Except the finest layer $$k=1$$, rest of $$a$$'s are real numbers.