[Python] Python preserve original matrix proportionality, scale with prescribed row sums and...

I have a historical sales mix profile (mxn) (rows = stores, columns = departments). I want to preserve this pattern/ profile as much as possible but set sales targets for every store and department (ie. row sums and columns sums predefined).

Method 1: I tried to scale rows and columns alternatively iteratively. But the targets don't seem to converge.

Method 2: I tried optimization algorithm scipy.optimize to minimize difference in matrix values proportions and constraint row sums and column sums but the it throws an error.

Any idea how I should correct my code?

Ps. non-negative matrix. If original value is 0, must preserve the 0.

Method 1:

def balance_matrix(original_matrix, target_store_sums, target_dept_sums, max_iterations=20, tolerance=1e-3):
"""
Adjust matrix values to match target store and dept sums while staying close to original values
and preserving zeros from the original matrix.

Args:
original_matrix: Initial 2D numpy array
target_store_sums: Array of target store sums
target_dept_sums: Array of target dept sums
max_iterations: Maximum number of iterations
tolerance: Convergence tolerance

Returns:
Balanced matrix meeting the constraints
"""

# Convert inputs to numpy arrays if they aren't already
matrix = np.array(original_matrix, dtype=float)
store_sums_target = np.array(target_store_sums, dtype=float)
dept_sums_target = np.array(target_dept_sums, dtype=float)

# Create mask for non-zero elements
mask = (original_matrix != 0)

# Verify dimensions
n_stores, n_depts = matrix.shape
assert len(store_sums_target) == n_stores
assert len(dept_sums_target) == n_depts

# Iterative scaling
for iteration in range(max_iterations):
prev_matrix = matrix.copy()

# Dept scaling
dept_sums = matrix.sum(axis=0)
dept_factors = np.ones(n_depts)
for j in range(n_depts):
if dept_sums[j] > 0:
# Only consider non-zero elements in scaling
dept_factors[j] = target_dept_sums[j] / dept_sums[j]

# Apply dept scaling only to non-zero elements
matrix = matrix * mask * dept_factors + matrix * (1 - mask)
matrix = np.maximum(matrix, 0)

# Store scaling
store_sums = matrix.sum(axis=1)
store_factors = np.ones(n_stores)
for i in range(n_stores):
if store_sums > 0:
# Only consider non-zero elements in scaling
store_factors = target_store_sums / store_sums

# Apply store scaling only to non-zero elements
matrix = matrix * mask * store_factors[:, np.newaxis] + matrix * (1 - mask)
matrix = np.maximum(matrix, 0)


#Check convergence
max_diff = np.max(np.abs(matrix - prev_matrix))
if max_diff < tolerance:
break

return matrix
matrix = np.array([[0, 2, 1], [3, 7, 8], [0, 1, 2]])
row_targets = [3, 3, 3]
col_targets = [0, 3, 3]
print(balance_matrix(matrix, row_targets, col_targets))


Method 2:

import numpy as np
from scipy.optimize import minimize
def optimize_matrix_quadratic(historical, row_targets, col_targets):
n_rows, n_cols = historical.shape

# Calculate initial proportions
total_sum = historical.sum()
proportions = historical / total_sum

# Flatten the matrix for optimization
x0 = historical.flatten()

def objective(x):
# Reshape x into matrix form
matrix = x.reshape((n_rows, n_cols))

# Calculate current proportions
current_proportions = matrix / matrix.sum()

w1, w2, w3 = 1.0, 3.0, 3.0 # Adjustable weights
# Penalty for deviating from original proportions
prop_penalty = w1*np.sum((current_proportions - proportions) ** 2)


return prop_penalty

# Constraints
constraints = [
{'type': 'eq', 'fun': np.sum(x, axis=1) - row_targets}, # Total sum constraint
{'type': 'eq', 'fun': np.sum(x, axis=0) - col_targets}, # Total sum constraint
]

# Non-negativity constraints
bounds = [(0, None) for _ in range(n_rows * n_cols)]

# Optimize
result = minimize(
objective,
x0,
method='SLSQP',
bounds=bounds,
constraints = constraints,
options={'maxiter': 1000}
)

# Reshape result back to matrix form
optimized = result.x.reshape((n_rows, n_cols))

return optimized


# Test the new methods
# Generate test data
np.random.seed(82)
historical = np.random.randint(2000, 15000, size=(6, 6))
row_targets = np.random.randint(50000, 70000, size=6)
col_targets = np.random.randint(45000, 65000, size=6)

# Run optimization using the hybrid approach
result = optimize_matrix_quadratic(historical, row_targets, col_targets)

print("Historical Matrix:")
print(historical)
print("\nOptimized Matrix (Hybrid):")
print(np.round(result).astype(int))

print("\nVerification:")
print("Row sums:", np.sum(result, axis=1).astype(int))
print("Row targets:", row_targets)
print("Column sums:", np.sum(result, axis=0).astype(int))
print("Column targets:", col_targets)

Continue reading...