Minimize sum of squares of rows in matrix when sum of columns have some constraint

Question

I'm looking for an algorithm that can find any matrix $a_{j,i}$ such that $$ \sum_{i \in I} \left(\sum_{j\in J} a_{j,i}\right)^2 $$ is minimal, while also for each $j\in J$ satisfying the constraint $$ k_j = \sum_{i \in I} a_{j,i} $$ where $k_j$ is some non-negative constant. Furthermore every $a_{j,i}$ must be non-negative, and some arbitrary selection of $a_{j,i}$s in the matrix are fixed to be zero.

---

This problem models a situation where each $j$ has some fixed amount of resources that need to be divided in some way between a row of buckets. The constraint saying that some values are fixed to zero means that some buckets are missing on some rows.

We then want to find the way to divide the resources between buckets such that the columns of buckets have their contents divided as equally as possible between them.

---

Notice that if no values are fixed to zero, the solution is easy, since we can simply assign the values $$ a_{j,i} = \frac{k_j}{\#I}, $$ where $\#I$ is the number of columns. This way every line of buckets has an equal amount of resources assigned, and since decreasing the contents in one line would increase the contents of another, this is optimal.

---

The only solution I've currently been able to find is using some sort of linear programming solver that supports quadratic objectives, but I hope that this problem could be solved more efficiently than that.

---

Nontrivial example Suppose we are looking for a 3x3 matrix where \begin{align} 1 &= a_{1,1} + a_{1,2} \\ 2 &= a_{2,1} + a_{2,2} + a_{2,3} \\ 1 &= \phantom{a_{3,1} +\,} a_{3,2} + a_{3,3} \end{align} Notice that some terms are missing in these sums. This is equivalent of fixing those terms to zero.

We then want to minimize $$ (a_{1,1}+a_{2,1})^2 + (a_{1,2}+a_{2,2}+a_{3,2})^2 + (a_{2,3}+a_{3,3})^2. $$ Notice that a lower bound on the minimum is $3\cdot\left(\frac{4}{3}\right)^2=\frac{16}4$, since in this case the $4$ total resources are equally divided.

In this case we can achieve that lower bound with the solution \begin{bmatrix} 1/3 & 2/3 & 0 \\ 1 & 0 & 1 \\ 0 & 2/3 & 1/3 \end{bmatrix} which indeed has the minimum $\frac{16}3$.

A different solution is \begin{bmatrix} 2/3 & 1/3 & 0 \\ 2/3 & 2/3 & 2/3 \\ 0 & 1/3 & 2/3 \end{bmatrix}

---

Example that doesn't hit lower bound of minimum. Suppose we have the constraints \begin{align} 1 &= a_{1,1} + a_{1,2} \\ 4 &= \phantom{a_{2,1} +}\,\, a_{2,2} \\ 1 &= \phantom{a_{3,1} +}\,\, a_{3,2} + a_{3,3} \end{align} In this case we have $a_{1,1} \leq 1$, but in order to reach the lower bound where every line is equal to $2$, we must have $a_{1,1} = 2$, which is not possible.

In this case there is only one optimal solution, which is: \begin{bmatrix} 1 & 0 & 0 \\ 0 & 4 & 0 \\ 0 & 0 & 1 \end{bmatrix}

Answer 1

Summary

Solve iteratively by pseudo-inverse, set negative $x_{ij}$ to 0, and remove negative $x_{ij}$ from equations.

Warning: seems to always work only if number of rows equal to the number of columns. Don't know if it is because of ill-conditioning of matrix or the initial guess is too far off or bug.

Lagrange Multipliers

I will write in math's ordering of index.

$x_{ij}$ is the element of matrix $X$ at row $i$ and column $j$.

The square of column sum of the $j^{th}$ column is:

$$\left(\sum_{i} x_{ij} \right)^{2}$$

The row sum constraint for the $i^{th}$ row is:

$$0 = s_i - \sum_{j} x_{ij} $$

Optimize sum of square of column sums by Lagrange multipliers:

$$ L = \left [ \sum_{j} \left(\sum_{i} x_{ij}\right)^2 \right ] - \sum_{i} \lambda_i\left(s_i - \sum_{j} x_{ij}\right) $$

Set derivatives to zero:

$$\forall i, j: 0 = \frac{\partial L}{\partial x_{ij}} = 2\left (\sum_{k} x_{kj}\right) + \lambda_{i}$$

$$\forall i: 0 = \frac{\partial L}{\partial \lambda_{i}} = s_i - \sum_{j}x_{ij} $$

Pseudo-inverse

Matrix equation:

$$A\vec{y} = \vec{b}$$

"Solve" by Moore-Penrose inverse (pseudo-inverse)

$$\vec{y} = A^{+}\vec{b}$$

This gives a solution that is as good as possible, even when some of the equations conflicts with each other.

Karush–Kuhn–Tucker conditions

What if the solution says some of the $x_{ij}$ is negative?

We set all negative $x_{ij}$ to zero, remove these variables from the matrix equation, and solve by Moore-Penrose inverse again. Repeat until no more negative $x_{ij}$.

(This is a simplified version of KKT because I don't understand everything in KKT yet.)

See: https://math.stackexchange.com/questions/1747435/correct-formulation-of-equality-and-non-negativity-constrained-non-linear-minimi

https://optimization.mccormick.northwestern.edu/index.php/Quadratic_programming

Example

I will explain with the 1-4-1 problem as an example.

Let $\vec{y}$ be the coefficients in row-major ordering and the $\lambda_i$

$$\vec{b} = \begin{bmatrix} x_{1,1}&x_{1,2}&x_{1,3}&x_{2,1}&x_{2,2}&x_{2,3}&x_{3,1}&x_{3,2}&x_{3,3}&\lambda_{1}/2&\lambda_{2}/2&\lambda_{3}/2\\ \end{bmatrix}^{T}$$

where $[\cdot]^{T}$ denotes transpose of a matrix/vector.

I define a geometry matrix $G$. If $g_{ij} = 0$, $x_{ij}$ should be zero. $$G = \begin{bmatrix} {\color{red}1}&{\color{red}1}&{\color{red}0}\\ {\color{green}0}&{\color{green}1}&{\color{green}0}\\ {\color{blue}0}&{\color{blue}1}&{\color{blue}1}\\ \end{bmatrix}$$

Matrix $A$ is:

$$A = \begin{bmatrix} {\color{red}1}&{\color{red}1}&{\color{red}0}&0&0&0&0&0&0&0&0&0\\ 0&0&0&{\color{green}0}&{\color{green}1}&{\color{green}0}&0&0&0&0&0&0\\ 0&0&0&0&0&0&{\color{blue}0}&{\color{blue}1}&{\color{blue}1}&0&0&0\\ \hdashline {\color{red}1}&0&0&{\color{green}0}&0&0&{\color{blue}0}&0&0&{\color{magenta}1}&0&0\\ 0&{\color{red}1}&0&0&{\color{green}1}&0&0&{\color{blue}1}&0&0&{\color{magenta}1}&0\\ 0&0&{\color{red}0}&0&0&{\color{green}0}&0&0&{\color{blue}1}&0&0&{\color{magenta}1}\\ \end{bmatrix}$$

• Upper half is for constraints.
• Bottom half is for $0 = {\partial L}/{\partial x_{ij}}$.
• The magenta numbers are for $\lambda_1$, $\lambda_2$, and $\lambda_3$.

The $\vec{b}$ is just the row sum padded by zeros:

$$\vec{b} = \begin{bmatrix}1&4&1&0&0&0\\\end{bmatrix}^{T}$$

Result

--- cycle 0 ---
    Matrix A:
        [[1. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
         [0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]
         [0. 0. 0. 0. 0. 0. 0. 1. 1. 0. 0. 0.]
         [1. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0.]
         [0. 1. 0. 0. 1. 0. 0. 1. 0. 0. 1. 0.]
         [0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 1.]]
    Vector b^{T}:
        [[1. 4. 1. 0. 0. 0.]]
    Vector y^{T}:
        [[ 1.4 -0.4  0.   0.   4.   0.   0.  -0.4  1.4 -1.4 -3.2 -1.4]]
    X:
        [[ 1.4 -0.4  0. ]
         [ 0.   4.   0. ]
         [ 0.  -0.4  1.4]]
    lambda^{T} / 2:
        [[-1.4 -3.2 -1.4]]
    row sum of X:
        [1. 4. 1.]
    column sum of X:
        [1.4 3.2 1.4]
--- cycle 1 ---
    Matrix A:
        [[1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
         [0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]
         [0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0.]
         [1. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0.]
         [0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 1. 0.]
         [0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 1.]]
    Vector b^{T}:
        [[1. 4. 1. 0. 0. 0.]]
    Vector y^{T}:
        [[ 1.  0.  0.  0.  4.  0.  0.  0.  1. -1. -4. -1.]]
    X:
        [[1. 0. 0.]
         [0. 4. 0.]
         [0. 0. 1.]]
    lambda^{T} / 2:
        [[-1. -4. -1.]]
    row sum of X:
        [1. 4. 1.]
    column sum of X:
        [1. 4. 1.]

Program

import numpy as np
from numpy.linalg import pinv, lstsq


def vec_y_split(vec_y, geometry):
    """Reshape vec_y into reshaped vec_y matrix and lambda
    """
    n_rows, n_cols = geometry.shape
    return vec_y[:-n_cols].reshape(geometry.shape), vec_y[-n_cols:]


def array_to_indent_str(x, indent=1):
    s = " " * indent * 4
    return s + str(x).replace("\n", "\n" + s)


def print_vec_y(vec_y, geometry, row_sum, indent=0):
    c, lam = vec_y_split(vec_y, geometry)
    s = " " * indent * 4
    print(s + "X:")
    print(array_to_indent_str(np.around(c, 3), indent + 1))
    print(s + "lambda^{T} / 2:")
    print(array_to_indent_str(np.around(lam.T, 3), indent + 1))
    print(s + "row sum of X:")
    print(array_to_indent_str(np.around(c.sum(axis=1), 3), indent + 1))
    print(s + "column sum of X:")
    print(array_to_indent_str(np.around(c.sum(axis=0), 3), indent + 1))
    print(s + "Error (Not considering zero requirement):")
    err = np.sum((c.sum(axis=1) - row_sum) ** 2) ** 0.5
    print(s + str(err))


def make_mat_A(geometry):
    n_rows, n_cols = geometry.shape
    mat_A = np.zeros((n_cols + n_rows, geometry.size + n_cols))
    # row equations
    for i in range(n_rows):
        for j in range(n_cols):
            mat_A[i, i * n_cols + j] = geometry[i, j]
    # col equations
    for j in range(n_cols):
        mat_A[n_rows + j, geometry.size + j] = 1
        for i in range(n_rows):
            mat_A[n_rows + j, i * n_cols + j] = geometry[i, j]
    return mat_A


def make_vec_b(geometry, row_sum):
    n_rows, n_cols = geometry.shape
    vec_b = np.empty((n_cols + n_rows, 1))
    vec_b[:n_rows, 0] = row_sum
    return vec_b


def remove_var(geometry, mat_A, coeff, epsilon):
    n_rows, n_cols = geometry.shape
    for i in range(n_rows):
        for j in range(n_cols):
            ix = i * n_cols + j
            if coeff[ix] < -epsilon:
                mat_A[:, ix] = 0


def iterate(geometry, mat_A, vec_b, vec_y, row_sum):
    i = 0
    epsilon = 1e-7
    n_rows, n_cols = geometry.shape

    while i == 0 or np.any(vec_y[:-n_cols] < -epsilon):
        print("--- cycle %d ---" % i)
        print("    Matrix A:")
        print(array_to_indent_str(np.around(mat_A, 3), 2))
        print("    Vector b^{T}:")
        print(array_to_indent_str(np.around(vec_b.T, 3), 2))

        vec_y = lstsq(mat_A, vec_b, rcond="warn")[0]

        if i > 0:
            remove_var(geometry, mat_A, vec_y, epsilon)

        print("    Vector y^{T}:")
        print(array_to_indent_str(np.around(vec_y.T, 3), 2))
        print_vec_y(vec_y, geometry, row_sum, indent=1)

        i += 1


def search(geometry, row_sum):
    print("=== Start ===")
    print("row_sum^{T}:")
    print(row_sum.T)
    print("geometry:")
    print(geometry)
    n_rows, n_cols = geometry.shape
    vec_y = np.empty(geometry.size + n_cols)
    # Set initial guess to solution with no zero
    vec_y[:-n_cols].reshape(n_rows, n_cols)[:] = \
        row_sum[:, np.newaxis] / n_cols
    vec_y[-n_cols:] = 0
    print("Intial Guess:")
    print_vec_y(vec_y, geometry, row_sum, indent=1)

    mat_A = make_mat_A(geometry)
    vec_b = make_vec_b(geometry, row_sum)
    iterate(geometry, mat_A, vec_b, vec_y, row_sum)


def problem_1_4_1():
    # For the 1-4-1 problem, use
    row_sum = np.array([1, 4, 1])
    geometry = np.array(
        [[1, 1, 0],
         [0, 1, 0],
         [0, 1, 1], ]
    )
    search(geometry, row_sum)


def problem_random():
    np.set_printoptions(precision=3)
    N = 100 * 2
    n_rows = N
    n_cols = N
    # row_sum is integers between 1 and 100 inclusive.
    row_sum = np.random.randint(1, 1000, size=n_rows)
    # geometry is a random array of zeros and ones.
    geometry = np.random.choice(2, size=(n_rows, n_cols))

    search(geometry, row_sum)


problem_1_4_1()
problem_random()

Random Implementation Ideas

1. No need to calculate pseudo inverse. Some solvers solves for a linear least square system. That is equivalent to solve by Moore-Penrose inverse, but possibly faster. https://eigen.tuxfamily.org/dox/group__LeastSquares.html

2. Use sparse solvers if there are many zeros.

3. Maybe rewrite the equations in block triangular form. For each overall iteration, perform one newton iteration on each block. That means no need for sparse matrix and two much smaller dense matrices. Don't know if that is possible though.

4. Use old solution as initial value for finding new solution.

Alternative and Questions

Use/Derive Simplex/criss-cross algorithm for the derivative of the lagrange multiplier equations. Don't know speed.

Or Use/Derive Simplex algorithm for quadratic programming.

https://www.me.utexas.edu/~jensen/ORMM/supplements/methods/nlpmethod/S2_quadratic.pdf

https://math.stackexchange.com/questions/246808/analog-of-simplex-method-for-quadratic-programming

What algorithm does the solver use?

How big is the matrix $X$?

How many zeros in matrix $X$?