def METHOD_NAME(
    self: Float[LinearOperator, "... N N"],
    right_tensor: Union[Float[Tensor, "... N P"], Float[Tensor, " N"]],
    left_tensor: Optional[Float[Tensor, "... O N"]] = None,
) -> Union[Float[Tensor, "... N P"], Float[Tensor, "... N"], Float[Tensor, "... O P"], Float[Tensor, "... O"]]:
    is_vector = right_tensor.ndim == 1
    if is_vector:
        right_tensor = right_tensor.unsqueeze(-1)
    res = self.root._cholesky_solve(right_tensor, upper=self.upper)
    if is_vector:
        res = res.squeeze(-1)
    if left_tensor is not None:
        res = left_tensor @ res
    return res