[BUG] probs doesn't work with Pauli words

[CatalinaAlbornoz]

Expected behavior

qml.probs() outputs the correct probabilities when using Pauli words as observables.

This is related to issue #6878 but in this case, because we're working with Pauli words, we may be able to find an easier solution.

Actual behavior

The probabilities are sometimes in the wrong order and sometimes just fully incorrect.

Additional information

For Pauli words, maybe we shouldn't need to take the tensor product of component eigvals, since we know a priori the correct eigenvalues and ordering of the full word.
E.g., this function https://docs.pennylane.ai/en/stable/code/api/pennylane.pauli.pauli_eigs.html should generate the correct eigenvalues in the correct order, the only input needed is the number of Pauli ops.
Optionally, you may then need to permute the list of eigenvalues if the wire order is not the same as the underlying device.

Source code

import pennylane as qml
import numpy as np

dev = qml.device("default.qubit", wires=3)

ob1 = qml.Y(0) @ qml.Y(1) @ qml.X(2)
ob2 = qml.X(2) @ qml.Y(0) @ qml.Y(1)

@qml.qnode(dev)
def circuit():
    qml.RX(0.4,wires=1)
    return qml.probs(op=ob1), qml.probs(op=ob2)

ans = circuit()
# In this case all the answers are different where they shouldn't be
print("ob1: ",ans[0],"ob2: ",ans[1],sep="\n")

Tracebacks

ob1: 
[0.07632271 0.07632271 0.17367729 0.17367729 0.07632271 0.07632271
 0.17367729 0.17367729]
ob2: 
[0.07632271 0.17367729 0.07632271 0.17367729 0.07632271 0.17367729
 0.07632271 0.17367729]

System information

Name: PennyLane
Version: 0.40.0
Summary: PennyLane is a cross-platform Python library for quantum computing, quantum machine learning, and quantum chemistry. Train a quantum computer the same way as a neural network.
Home-page: https://github.com/PennyLaneAI/pennylane
Author: 
Author-email: 
License: Apache License 2.0
Location: /usr/local/lib/python3.11/dist-packages
Requires: appdirs, autograd, autoray, cachetools, diastatic-malt, networkx, numpy, packaging, pennylane-lightning, requests, rustworkx, scipy, tomlkit, typing-extensions
Required-by: PennyLane_Lightning

Platform info:           Linux-6.1.85+-x86_64-with-glibc2.35
Python version:          3.11.11
Numpy version:           1.26.4
Scipy version:           1.13.1
Installed devices:
- default.clifford (PennyLane-0.40.0)
- default.gaussian (PennyLane-0.40.0)
- default.mixed (PennyLane-0.40.0)
- default.qubit (PennyLane-0.40.0)
- default.qutrit (PennyLane-0.40.0)
- default.qutrit.mixed (PennyLane-0.40.0)
- default.tensor (PennyLane-0.40.0)
- null.qubit (PennyLane-0.40.0)
- reference.qubit (PennyLane-0.40.0)
- lightning.qubit (PennyLane_Lightning-0.40.0)

Existing GitHub issues

• I have searched existing GitHub issues to make sure the issue does not already exist.

[CatalinaAlbornoz]

Analysis by @JerryChen97 :

It looks like there's fundamental obstacle to proceed:

Essentially, the qml.probs interface we provided when an op= was given, is not actually a 'physical' output that is invariant. Instead, it gives the probability of the rotated state, rotation matrix of which was given by the diagonalization of the given operator. For example, if op=Y , then in PennyLane, we find the U such that U^dagger Y U = -Z, which turns out to be U=Z S H , and then apply the U to the initial state to get the rotated state and finally provide the probs of that state.

However, it is well known that the eigen-decomposition of an Hermitian matrix has a huge gauge freedom. Still the same example, one can choose U=S H which made life much easier since now Y -> Z. For operators, this gauge doesn't matter at all because whatever decomposition was given, it is nothing different but in another basis set. But this matters when it comes to the state that gets applied in PL pipeline, because a state only gets applied ONCE by the U! Different U brings completely different state_applied_by_U , especially when the original operator has a degenerated eigen-subspace (e.g. X@Y, having deg=2 for both +1 and -1), there will be multiple independent U(n) freedom all tensor producted together.
When wires=1 this can still be walked around by replacing the eigh in Hermitian.eigendecomposition with eig, which will keep the eigenvalues' order the same as qml.H's intrinsic decomposition no sorting atl. BUT when it comes to multiple wires, this is almost impossible to resolve because the freedom is too much huge and how can we guarantee we always have the same 'standard' order?...

Suggestion:
Add a UserWarning to probs when non-trivial op option used.
For now only Z and I are considered trivial.
Alternatively, we can first unify the eigvals orders of all the single-qubit operators, then consider them as trivial as well.