Minimize sum of squares of rows in matrix when sum of columns have some constraint
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}&amp;x_{1,2}&amp;x_{1,3}&amp;x_{2,1}&amp;x_{2,2}&amp;x_{2,3}&amp;x_{3,1}&amp;x_{3,2}&amp;x_{3,3}&amp;\lambda_{1}/2&amp;\lambda_{2}/2&amp;\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}&amp;{\color{red}1}&amp;{\color{red}0}\\
{\color{green}0}&amp;{\color{green}1}&amp;{\color{green}0}\\
{\color{blue}0}&amp;{\color{blue}1}&amp;{\color{blue}1}\\
\end{bmatrix}$$

Matrix $A$ is:

$$A = \begin{bmatrix}
{\color{red}1}&amp;{\color{red}1}&amp;{\color{red}0}&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0\\
0&amp;0&amp;0&amp;{\color{green}0}&amp;{\color{green}1}&amp;{\color{green}0}&amp;0&amp;0&amp;0&amp;0&amp;0&amp;0\\
0&amp;0&amp;0&amp;0&amp;0&amp;0&amp;{\color{blue}0}&amp;{\color{blue}1}&amp;{\color{blue}1}&amp;0&amp;0&amp;0\\ \hdashline
{\color{red}1}&amp;0&amp;0&amp;{\color{green}0}&amp;0&amp;0&amp;{\color{blue}0}&amp;0&amp;0&amp;{\color{magenta}1}&amp;0&amp;0\\
0&amp;{\color{red}1}&amp;0&amp;0&amp;{\color{green}1}&amp;0&amp;0&amp;{\color{blue}1}&amp;0&amp;0&amp;{\color{magenta}1}&amp;0\\
0&amp;0&amp;{\color{red}0}&amp;0&amp;0&amp;{\color{green}0}&amp;0&amp;0&amp;{\color{blue}1}&amp;0&amp;0&amp;{\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&amp;4&amp;1&amp;0&amp;0&amp;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] &lt; -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] &lt; -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 &gt; 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


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
Use sparse solvers if there are many zeros.
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.
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$?