[M-Colley]

Summary

Multiplying a batched MulLinearOperator by a batch of constants (a *batch x 1 x 1 tensor) fails, because LinearOperator.mul passes the *batch-shaped constants to _mul_constant, and MulLinearOperator._mul_constant tests them with if other > 0:

import torch
from linear_operator.operators import DenseLinearOperator, MulLinearOperator

a, b = torch.randn(2, 3, 4, 4).unbind(0)
A, B = a @ a.mT + torch.eye(4), b @ b.mT + torch.eye(4)
c = torch.tensor([2.0, 3.0, 0.5]).view(3, 1, 1)
MulLinearOperator(DenseLinearOperator(A), DenseLinearOperator(B)) * c
# RuntimeError: Boolean value of Tensor with more than one value is ambiguous

In BoTorch, this happens for a multi-output SingleTaskMultiFidelityGP with two data fidelities at a b x 1 x d batch of points:

import torch
from botorch.models import SingleTaskMultiFidelityGP

torch.set_default_dtype(torch.float64)
m = SingleTaskMultiFidelityGP(torch.rand(8, 4), torch.rand(8, 2), data_fidelities=[2, 3])
m.posterior(torch.rand(4, 1, 4))
# RuntimeError: Boolean value of Tensor with more than one value is ambiguous

With this PR, the posterior equals the posteriors computed point by point.

Found by the equivalence tests of botorch-rust.

Fix

• MulLinearOperator._mul_constant uses (other > 0).all(), as RootLinearOperator._mul_constant does. If all constants are positive, they are still pushed into the left root decomposition; otherwise the existing ConstantMulLinearOperator fallback is used.
• TriangularLinearOperator._mul_constant multiplied its *batch x N x N tensor by other.unsqueeze(-1), which scales the rows rather than the batch members. The result is wrong if the batch size equals N and an error otherwise. It is reached from the first change through RootLinearOperator._mul_constant when the roots are Cholesky factors (the default root decomposition of small dense matrices). It now uses other.unsqueeze(-1).unsqueeze(-1).

from linear_operator.operators import TriangularLinearOperator

L = torch.randn(4, 4, 4).tril()
c = torch.rand(4, 1, 1) + 0.5
torch.allclose((TriangularLinearOperator(L) * c).to_dense(), L * c)
# False before, True after

Scalar constants behave as before.

Not changed here: #148 makes the same change to LowRankRootAddedDiagLinearOperator._mul_constant, and #136 fixes batches of constants for the block operators. ConstantDiagLinearOperator._mul_constant multiplies its *batch x 1 values by the *batch constants without unsqueezing them, so AddedDiagLinearOperator(MulLinearOperator(...), ConstantDiagLinearOperator(...)) * c still fails for a batch of constants.

Tests

• test_constant_mul_batch in test/operators/test_mul_linear_operator.py (batch and multi-batch classes): positive and mixed-sign batches of constants, for MulLinearOperators of RootLinearOperators and of dense matrices with Cholesky roots. Positive constants keep the MulLinearOperator; the dense results match.
• test_constant_mul_batch in test/operators/test_triangular_linear_operator.py (lower and upper): batch shapes (5,) (equal to N) and (2, 3).

All four fail on main and pass with this PR.

python -m pytest test/operators/test_mul_linear_operator.py test/operators/test_triangular_linear_operator.py \
    test/operators/test_root_linear_operator.py test/operators/test_chol_linear_operator.py \
    test/operators/test_constant_mul_linear_operator.py test/operators/test_sum_linear_operator.py
# 1171 passed

ufmt check (black 22.3.0, usort 1.0.3) and flake8 pass on the changed files.

🤖 Generated with Claude Code