[Python] scipy 1.16 bug, invalid matrices with incorrect indptr sizes

I have just upgraded from 1.15 to 1.16, but rand across an issue with CSC matrices and multiplication.

Demonstration of scipy sparse matrix multiply bug affecting CSC matrix indptr size.

This script demonstrates a bug in scipy where the .multiply() operation on CSC matrices can produce invalid CSC matrices with incorrect indptr array sizes. This bug was first noticed in scipy 1.16.0 and affects downstream operations that rely on the CSC format invariant: len(indptr) == n_columns + 1.

The bug manifests when:

1. Two CSC matrices are multiplied using .multiply()
2. The result has len(indptr) != shape[1] + 1
3. This breaks downstream operations like matrix reconstruction

This script provides a minimal reproducible example and shows the fix.

"""
Demonstration of scipy sparse matrix multiply bug affecting CSC matrix indptr size.

This script demonstrates a bug in scipy where the .multiply() operation on CSC matrices
can produce invalid CSC matrices with incorrect indptr array sizes. This bug was first
noticed in scipy 1.16.0 and affects downstream operations that rely on the CSC format
invariant: len(indptr) == n_columns + 1.

The bug manifests when:
1. Two CSC matrices are multiplied using .multiply()
2. The result has len(indptr) != shape[1] + 1
3. This breaks downstream operations like matrix reconstruction

This script provides a minimal reproducible example and shows the fix.
"""

import numpy as np
from scipy.sparse import coo_matrix, csc_matrix
import scipy
import sys

def get_version_info():
"""Get and display version information."""
print("=" * 60)
print("SCIPY MULTIPLY BUG DEMONSTRATION")
print("=" * 60)
print(f"Python version: {sys.version}")
print(f"NumPy version: {np.__version__}")
print(f"SciPy version: {scipy.__version__}")
print("=" * 60)
print()

def create_test_matrices():
"""
Create test matrices that trigger the scipy multiply bug.

These matrices are designed to produce the bug when multiplied.
The data comes from actual sparse matrix operations that exhibited this issue.
"""
print("Creating test matrices...")

# Matrix 1 data (example sparse matrix)
# Raw data: [7, 0, 8, 2, 9, 3, 1, 1, 2, 2, 3, 3, 0, 1, 3, 2, 3, 4, 3, 1, 3, 0, 4, 4, 2, 4, 5, 4, 2, 3, 4, 3, 5, 0, 8, 4, 5, 0, 7, 0, 3, 4]
# Reshaped to (row_indices, col_indices) pairs
b0_raw = np.array([7, 0, 8, 2, 9, 3, 1, 1, 2, 2, 3, 3, 0, 1, 3, 2, 3, 4, 3, 1, 3, 0, 4, 4, 2, 4, 5, 4, 2, 3, 4, 3, 5, 0, 8, 4, 5, 0, 7, 0, 3, 4])
b0_reshaped = b0_raw.reshape((int(len(b0_raw) / 2), 2)).T

# Matrix 2 data (example sparse matrix)
# Raw data: [1, 1, 8, 3, 4, 2, 9, 3, 5, 1, 7, 4, 1, 1, 2, 0, 0, 0, 4, 4, 1, 3, 0, 3, 5, 1, 5, 0, 7, 4, 4, 0, 5, 3]
# Reshaped to (row_indices, col_indices) pairs
b1_raw = np.array([1, 1, 8, 3, 4, 2, 9, 3, 5, 1, 7, 4, 1, 1, 2, 0, 0, 0, 4, 4, 1, 3, 0, 3, 5, 1, 5, 0, 7, 4, 4, 0, 5, 3])
b1_reshaped = b1_raw.reshape((int(len(b1_raw) / 2), 2)).T

print(f"Matrix 1 shape: {b0_reshaped.shape}")
print(f"Matrix 2 shape: {b1_reshaped.shape}")
print()

# Create sparse matrices in CSC format
# Use ones for data (incidence matrices)
mat1 = coo_matrix((np.ones(b0_reshaped.shape[1]), b0_reshaped), shape=(10, 5)).tocsc()
mat2 = coo_matrix((np.ones(b1_reshaped.shape[1]), b1_reshaped), shape=(10, 5)).tocsc()

return mat1, mat2

def analyze_matrix(matrix, name):
"""
Analyze a sparse matrix and display its properties.

Args:
matrix: The sparse matrix to analyze
name: Name for display purposes
"""
print(f"{name}:")
print(f" Shape: {matrix.shape}")
print(f" nnz (non-zero elements): {matrix.nnz}")
print(f" indptr size: {len(matrix.indptr)}")
print(f" indices size: {len(matrix.indices)}")
print(f" data size: {len(matrix.data)}")
print(f" Expected indptr size: {matrix.shape[1] + 1}")
print(f" indptr valid: {len(matrix.indptr) == matrix.shape[1] + 1}")
print(f" indptr: {matrix.indptr}")
print(f" indices: {matrix.indices}")
print(f" data: {matrix.data}")
print()

def demonstrate_bug(mat1, mat2):
"""
Demonstrate the scipy multiply bug.

This function shows how .multiply() can produce an invalid CSC matrix
with incorrect indptr size in certain scipy versions.
"""
print("=" * 60)
print("DEMONSTRATING THE BUG")
print("=" * 60)

# Show original matrices
print("Original matrices:")
analyze_matrix(mat1, "Matrix 1")
analyze_matrix(mat2, "Matrix 2")

# Perform the problematic operation
print("Performing element-wise multiplication...")
result = mat1.multiply(mat2)

print("Result after .multiply():")
analyze_matrix(result, "Result")

# Check if the result is valid
is_valid = len(result.indptr) == result.shape[1] + 1
print(f"Result matrix is valid: {is_valid}")

if not is_valid:
print("❌ BUG DETECTED: Invalid indptr size!")
print(f" Expected: {result.shape[1] + 1}, Got: {len(result.indptr)}")
else:
print("✅ No bug detected: Valid indptr size")

print()

def test_matrix_reconstruction(matrix, name):
"""
Test if a matrix can be reconstructed from its components.

This simulates what happens when loading matrix data from storage.
"""
print(f"Testing reconstruction of {name}...")
try:
reconstructed = csc_matrix((matrix.data, matrix.indices, matrix.indptr),
shape=matrix.shape)
print(f" ✅ Reconstruction successful")
print(f" Reconstructed shape: {reconstructed.shape}")
print(f" Reconstructed nnz: {reconstructed.nnz}")
print(f" Reconstructed indptr size: {len(reconstructed.indptr)}")
except Exception as e:
print(f" ❌ Reconstruction failed: {e}")
print()

def demonstrate_fix(mat1, mat2):
"""
Demonstrate the fix for the scipy multiply bug.

The fix converts the invalid CSC matrix to COO format and back to CSC,
which reconstructs the matrix with correct indptr size.
"""
print("=" * 60)
print("DEMONSTRATING THE FIX")
print("=" * 60)

# Perform the problematic operation
result_broken = mat1.multiply(mat2)

print("Before fix:")
analyze_matrix(result_broken, "Broken Result")

# Apply the fix
if len(result_broken.indptr) != result_broken.shape[1] + 1:
print("Applying fix: .tocoo().tocsc()")
result_fixed = result_broken.tocoo().tocsc()
else:
print("No fix needed: Matrix is already valid")
result_fixed = result_broken

print("After fix:")
analyze_matrix(result_fixed, "Fixed Result")

# Test reconstruction
test_matrix_reconstruction(result_fixed, "Fixed Matrix")

# Show dense matrix for comparison
print("Dense matrix representation:")
dense_result = result_fixed.toarray()
print(dense_result)
print()

def show_dense_comparison(mat1, mat2):
"""
Show the dense matrix comparison to prove the mathematical result is correct.
"""
print("=" * 60)
print("DENSE MATRIX COMPARISON")
print("=" * 60)

print("Matrix 1 (dense):")
print(mat1.toarray())
print()

print("Matrix 2 (dense):")
print(mat2.toarray())
print()

# Element-wise multiplication in dense format
dense_result = mat1.toarray() * mat2.toarray()
print("Element-wise multiplication result (dense):")
print(dense_result)
print()

# Show sparse result for comparison
sparse_result = mat1.multiply(mat2)
print("Element-wise multiplication result (sparse):")
print(sparse_result.toarray())
print()

# Check if they match
dense_sparse_match = np.allclose(dense_result, sparse_result.toarray())
print(f"Dense and sparse results match: {dense_sparse_match}")
print()

def main():
"""
Main function to run the complete demonstration.
"""
get_version_info()

# Create test matrices
mat1, mat2 = create_test_matrices()

# Demonstrate the bug
demonstrate_bug(mat1, mat2)

# Show dense comparison
show_dense_comparison(mat1, mat2)

# Demonstrate the fix
demonstrate_fix(mat1, mat2)

print("=" * 60)
print("SUMMARY")
print("=" * 60)
print("This script demonstrates a bug in scipy where .multiply() on CSC matrices")
print("can produce invalid matrices with incorrect indptr sizes.")
print()
print("The bug affects:")
print("- Matrix reconstruction from components")
print("- Matrix storage and loading")
print("- Any operation that relies on CSC format invariants")
print()
print("The fix (.tocoo().tocsc()) restores the matrix to valid CSC format")
print("while preserving the mathematical result.")
print()
print("This bug was first noticed in scipy 1.16.0 and may affect other versions.")
print("=" * 60)

if __name__ == "__main__":
main()


And here is the output:

============================================================
SCIPY MULTIPLY BUG DEMONSTRATION
============================================================
Python version: 3.12.0 (main, Mar 15 2024, 15:20:13) [Clang 15.0.0 (clang-1500.1.0.2.5)]
NumPy version: 1.26.4
SciPy version: 1.16.0
============================================================

Creating test matrices...
Matrix 1 shape: (2, 21)
Matrix 2 shape: (2, 17)

============================================================
DEMONSTRATING THE BUG
============================================================
Original matrices:
Matrix 1:
Shape: (10, 5)
nnz (non-zero elements): 18
indptr size: 6
indices size: 18
data size: 18
Expected indptr size: 6
indptr valid: True
indptr: [ 0 3 6 9 13 18]
indices: [3 5 7 0 1 3 2 3 8 2 3 4 9 2 3 4 5 8]
data: [1. 2. 2. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 2. 1. 1. 1.]

Matrix 2:
Shape: (10, 5)
nnz (non-zero elements): 14
indptr size: 6
indices size: 14
data size: 14
Expected indptr size: 6
indptr valid: True
indptr: [ 0 4 6 7 12 14]
indices: [0 2 4 5 1 5 4 0 1 5 8 9 4 7]
data: [1. 1. 1. 1. 2. 2. 1. 1. 1. 1. 1. 1. 1. 2.]

Performing element-wise multiplication...
Result after .multiply():
Result:
Shape: (10, 5)
nnz (non-zero elements): 4
indptr size: 11
indices size: 4
data size: 4
Expected indptr size: 6
indptr valid: False
indptr: [0 0 1 1 1 2 3 3 3 3 4]
indices: [1 4 0 3]
data: [2. 1. 2. 1.]

Result matrix is valid: False
❌ BUG DETECTED: Invalid indptr size!
Expected: 6, Got: 11

============================================================
DENSE MATRIX COMPARISON
============================================================
Matrix 1 (dense):
[[0. 1. 0. 0. 0.]
[0. 1. 0. 0. 0.]
[0. 0. 1. 1. 1.]
[1. 1. 1. 1. 2.]
[0. 0. 0. 1. 1.]
[2. 0. 0. 0. 1.]
[0. 0. 0. 0. 0.]
[2. 0. 0. 0. 0.]
[0. 0. 1. 0. 1.]
[0. 0. 0. 1. 0.]]

Matrix 2 (dense):
[[1. 0. 0. 1. 0.]
[0. 2. 0. 1. 0.]
[1. 0. 0. 0. 0.]
[0. 0. 0. 0. 0.]
[1. 0. 1. 0. 1.]
[1. 2. 0. 1. 0.]
[0. 0. 0. 0. 0.]
[0. 0. 0. 0. 2.]
[0. 0. 0. 1. 0.]
[0. 0. 0. 1. 0.]]

Element-wise multiplication result (dense):
[[0. 0. 0. 0. 0.]
[0. 2. 0. 0. 0.]
[0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0.]
[0. 0. 0. 0. 1.]
[2. 0. 0. 0. 0.]
[0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0.]
[0. 0. 0. 1. 0.]]

Element-wise multiplication result (sparse):
[[0. 0. 0. 0. 0.]
[0. 2. 0. 0. 0.]
[0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0.]
[0. 0. 0. 0. 1.]
[2. 0. 0. 0. 0.]
[0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0.]
[0. 0. 0. 1. 0.]]

Dense and sparse results match: True

============================================================
DEMONSTRATING THE FIX
============================================================
Before fix:
Broken Result:
Shape: (10, 5)
nnz (non-zero elements): 4
indptr size: 11
indices size: 4
data size: 4
Expected indptr size: 6
indptr valid: False
indptr: [0 0 1 1 1 2 3 3 3 3 4]
indices: [1 4 0 3]
data: [2. 1. 2. 1.]

Applying fix: .tocoo().tocsc()
After fix:
Fixed Result:
Shape: (10, 5)
nnz (non-zero elements): 4
indptr size: 6
indices size: 4
data size: 4
Expected indptr size: 6
indptr valid: True
indptr: [0 1 2 2 3 4]
indices: [5 1 9 4]
data: [2. 2. 1. 1.]

Testing reconstruction of Fixed Matrix...
✅ Reconstruction successful
Reconstructed shape: (10, 5)
Reconstructed nnz: 4
Reconstructed indptr size: 6

Dense matrix representation:
[[0. 0. 0. 0. 0.]
[0. 2. 0. 0. 0.]
[0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0.]
[0. 0. 0. 0. 1.]
[2. 0. 0. 0. 0.]
[0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0.]
[0. 0. 0. 1. 0.]]

============================================================
SUMMARY
============================================================
This script demonstrates a bug in scipy where .multiply() on CSC matrices
can produce invalid matrices with incorrect indptr sizes.

The bug affects:
- Matrix reconstruction from components
- Matrix storage and loading
- Any operation that relies on CSC format invariants

The fix (.tocoo().tocsc()) restores the matrix to valid CSC format
while preserving the mathematical result.

This bug was first noticed in scipy 1.16.0 and may affect other versions.
============================================================


Prior versions work fine.

Continue reading...