Repository: PennyLaneAI/pennylane
Pull request #6776 (2025-01-07 17:11:14)

matrix norm method
norm computations

github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [`doc/releases/changelog-dev.md`](/PennyLaneAI/pennylane/blob/master/doc/releases/changelog-dev.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.

Patch:
diff --git a/pennylane/labs/vibronic/v2/matrix_ops.py b/pennylane/labs/vibronic/v2/matrix_ops.py
index a3c87aeab..c0746ccd4 100644
--- a/pennylane/labs/vibronic/v2/matrix_ops.py
+++ b/pennylane/labs/vibronic/v2/matrix_ops.py
@@ -136,3 +136,8 @@ def _zeros(shape: Tuple[int], sparse: bool = False):
         return sp.sparse.csr_matrix(shape)
 
     return np.zeros(shape, dtype=np.complex128)
+
+
+def op_norm(gridpoints: int) -> float:
+    """The norm of P and Q"""
+    return np.sqrt(gridpoints * np.pi / 2)
diff --git a/pennylane/labs/vibronic/v2/utils.py b/pennylane/labs/vibronic/v2/utils.py
index cafee906a..eba5c021c 100644
--- a/pennylane/labs/vibronic/v2/utils.py
+++ b/pennylane/labs/vibronic/v2/utils.py
@@ -3,7 +3,7 @@
 import numpy as np
 
 
-def coeffs(states, modes):
+def coeffs(states: int, modes: int):
     """Produce random coefficients for input"""
 
     alphas = np.random.random(size=(states, states, modes))
@@ -26,3 +26,11 @@ def coeffs(states, modes):
             symmetric_lambdas[i, j] = (lambdas[i, j] + lambdas[j, i]) / 2
 
     return symmetric_alphas, symmetric_betas, symmetric_lambdas, omegas
+
+
+def is_pow_2(k: int) -> bool:
+    """Test if k is a power of two"""
+    if (k & (k - 1) == 0) or k == 0:
+        return True
+
+    return False
diff --git a/pennylane/labs/vibronic/v2/vibronic_hamiltonian.py b/pennylane/labs/vibronic/v2/vibronic_hamiltonian.py
index 45efcf524..5c5977187 100644
--- a/pennylane/labs/vibronic/v2/vibronic_hamiltonian.py
+++ b/pennylane/labs/vibronic/v2/vibronic_hamiltonian.py
@@ -4,6 +4,7 @@ from __future__ import annotations
 
 import numpy as np
 from scipy.sparse import csr_matrix
+from utils import is_pow_2
 from vibronic_matrix import VibronicMatrix, commutator
 from vibronic_term import VibronicTerm, VibronicWord
 from vibronic_tree import Node
@@ -21,14 +22,23 @@ class VibronicHamiltonian:
         lambdas: np.ndarray,
         omegas: np.ndarray,
     ):
+        if not is_pow_2(states) or states == 0:
+            raise ValueError(f"States must be a positive power of 2, got {states} states.")
+
         if alphas.shape != (states, states, modes):
-            raise TypeError
+            raise TypeError(
+                f"Alphas must have shape {(states, states, modes)}, got shape {alphas.shape}."
+            )
         if betas.shape != (states, states, modes, modes):
-            raise TypeError
+            raise TypeError(
+                f"Betas must have shape {(states, states, modes, modes)}, got shape {betas.shape}."
+            )
         if lambdas.shape != (states, states):
-            raise TypeError
+            raise TypeError(
+                f"Lambdas must have shape {(states, states)}, got shape {lambdas.shape}."
+            )
         if omegas.shape != (modes,):
-            raise TypeError
+            raise TypeError(f"Omegas must have shape {(modes,)}, got shape {omegas.shape}.")
 
         self.states = states
         self.modes = modes
diff --git a/pennylane/labs/vibronic/v2/vibronic_matrix.py b/pennylane/labs/vibronic/v2/vibronic_matrix.py
index 936af9fbd..9678c3a82 100644
--- a/pennylane/labs/vibronic/v2/vibronic_matrix.py
+++ b/pennylane/labs/vibronic/v2/vibronic_matrix.py
@@ -2,12 +2,14 @@
 
 from __future__ import annotations
 
+import math
 from itertools import product
 from typing import Dict, Tuple, Union
 
 import numpy as np
 import scipy as sp
-from matrix_ops import _kron, _zeros, word_to_matrix
+from matrix_ops import _kron, _zeros, op_norm, word_to_matrix
+from utils import is_pow_2
 from vibronic_term import VibronicWord
 
 
@@ -17,6 +19,10 @@ class VibronicMatrix:
     def __init__(
         self, states: int, modes: int, blocks: Dict[Tuple[int, int], VibronicWord] = None
     ) -> VibronicMatrix:
+
+        if not is_pow_2(states) or states == 0:
+            raise ValueError(f"The number of states must be a power of 2, got {states}.")
+
         if blocks is None:
             blocks = {}
 
@@ -26,12 +32,24 @@ class VibronicMatrix:
 
     def block(self, row: int, col: int) -> VibronicWord:
         """Return the block indexed at (row, col)"""
+        if row < 0 or col < 0:
+            raise IndexError(f"Index cannot be negative, got {(row, col)}.")
+        if row >= self.states or col >= self.states:
+            raise IndexError(
+                f"Index out of bounds. Got {(row, col)} but there are only {self.states} states."
+            )
         return self._blocks.get((row, col), VibronicWord(tuple()))
 
     def set_block(self, row: int, col: int, word: VibronicWord) -> None:
         """Set the value of the block indexed at (row, col)"""
         if not isinstance(word, VibronicWord):
             raise TypeError(f"Block value must be VibronicWord. Got {type(word)}.")
+        if row < 0 or col < 0:
+            raise IndexError(f"Index cannot be negative, got {(row, col)}.")
+        if row >= self.states or col >= self.states:
+            raise IndexError(
+                f"Index out of bounds. Got {(row, col)} but there are only {self.states} states."
+            )
 
         self._blocks[(row, col)] = word
 
@@ -57,6 +75,52 @@ class VibronicMatrix:
 
         return matrix
 
+    def norm(self, gridpoints: int) -> float:
+        """Compute the spectral norm"""
+
+        if not is_pow_2(gridpoints) or gridpoints <= 0:
+            raise ValueError(
+                f"Number of gridpoints must be a positive power of 2, got {gridpoints}."
+            )
+
+        return self._norm(gridpoints)
+
+    def _norm(self, gridpoints: int) -> float:
+        # pylint: disable=protected-access
+        if len(self._blocks) == 0:
+            return 0
+
+        if self.states == 1:
+            return self._norm_base_case(gridpoints)
+
+        top_left, top_right, bottom_left, bottom_right = self._partition_into_quadrants()
+
+        norm1 = max(top_left._norm(gridpoints), bottom_right._norm(gridpoints))
+        norm2 = np.sqrt(top_right._norm(gridpoints) * bottom_left._norm(gridpoints))
+
+        return norm1 + norm2
+
+    def _norm_base_case(self, gridpoints: int) -> float:
+        if self.states != 1:
+            raise RuntimeError("Base case called on VibronicMatrix with >1 state")
+
+        if len(self._blocks) == 0:
+            return 0
+
+        word = self._blocks[(0, 0)]
+        norm = 0
+
+        for term in word.terms:
+            term_op_norm = math.prod(map(lambda op: op_norm(gridpoints) ** len(op), term.ops))
+            coeff_sum = sum(
+                term.coeffs.compute(index)
+                for index in product(range(self.modes), repeat=len(term.ops))
+            )
+
+            norm += coeff_sum * term_op_norm
+
+        return norm
+
     def __add__(self, other: VibronicMatrix) -> VibronicMatrix:
         if self.states != other.states:
             raise ValueError(
@@ -150,6 +214,29 @@ class VibronicMatrix:
 
         return True
 
+    def _partition_into_quadrants(self) -> Tuple[VibronicMatrix]:
+        # pylint: disable=chained-comparison
+        half = self.states // 2
+
+        top_left = VibronicMatrix(half, self.modes, {})
+        top_right = VibronicMatrix(half, self.modes, {})
+        bottom_left = VibronicMatrix(half, self.modes, {})
+        bottom_right = VibronicMatrix(half, self.modes, {})
+
+        for index, word in self._blocks.items():
+            x, y = index
+
+            if x < half and y < half:
+                top_left.set_block(x, y, word)
+            if x < half and y >= half:
+                top_right.set_block(x, y - half, word)
+            if x >= half and y < half:
+                bottom_left.set_block(x - half, y, word)
+            if x >= half and y >= half:
+                bottom_right.set_block(x - half, y - half, word)
+
+        return top_left, top_right, bottom_left, bottom_right
+
 
 def commutator(a: VibronicMatrix, b: VibronicMatrix):
     """Return the commutator [a, b] = ab - ba"""
diff --git a/pennylane/labs/vibronic/v2/vibronic_tree.py b/pennylane/labs/vibronic/v2/vibronic_tree.py
index dd1b65cc1..ad08d9f07 100644
--- a/pennylane/labs/vibronic/v2/vibronic_tree.py
+++ b/pennylane/labs/vibronic/v2/vibronic_tree.py
@@ -193,7 +193,7 @@ class Node:  # pylint: disable=too-many-instance-attributes
             ret += str(self.r_child, level + 1)
             ret += ")"
 
-        if self.node_type == NodeType.SUM:
+        if self.node_type == NodeType.HADAMARD:
             ret += "(HADAMARD, \n"
             ret += str(self.l_child, level + 1) + ",\n"
             ret += str(self.r_child, level + 1)
@@ -205,7 +205,7 @@ class Node:  # pylint: disable=too-many-instance-attributes
         """Compute the coefficient at the given index"""
 
         if not self._validate_index(index):
-            raise ValueError(f"Given index {index} is not of shape {self.shape}")
+            raise ValueError(f"Given index {index} is not compatible with shape {self.shape}")
 
         if self.node_type == NodeType.TENSOR:
             return self.tensor[index]


Test patch:
diff --git a/pennylane/labs/vibronic/v2/test_matrix.py b/pennylane/labs/vibronic/v2/test_matrix.py
new file mode 100644
index 000000000..0c64be4a3
--- /dev/null
+++ b/pennylane/labs/vibronic/v2/test_matrix.py
@@ -0,0 +1,70 @@
+"""Test matrix operations"""
+
+import numpy as np
+import pytest
+from matrix_ops import op_norm
+from vibronic_matrix import VibronicMatrix
+from vibronic_term import VibronicTerm, VibronicWord
+from vibronic_tree import Node
+
+vword0 = VibronicWord([VibronicTerm(tuple(), Node.tensor_node(np.array(0.5)))])
+blocks0 = {(0, 0): vword0}
+vmat0 = VibronicMatrix(1, 1, blocks0)
+
+vword1 = VibronicWord([VibronicTerm(("P",), Node.tensor_node(np.array([1])))])
+blocks1 = {(0, 0): vword1}
+vmat1 = VibronicMatrix(1, 1, blocks1)
+
+vword2a = VibronicWord([VibronicTerm(("P", "P"), Node.tensor_node(np.array([[0, 1], [2, 3]])))])
+vword2b = VibronicWord([VibronicTerm(("P",), Node.tensor_node(np.array([1, 2])))])
+blocks2 = {(0, 0): vword2a, (1, 1): vword2b}
+vmat2 = VibronicMatrix(2, 2, blocks2)
+
+blocks3 = {(0, 1): vword2a, (1, 0): vword2b}
+vmat3 = VibronicMatrix(2, 2, blocks3)
+
+blocks4 = {(0, 0): vword2a, (0, 1): vword2a, (1, 0): vword2b, (1, 1): vword2b}
+vmat4 = VibronicMatrix(2, 2, blocks4)
+
+
+class TestMatrix:
+    """Test properties of the VibronicMatrix class"""
+
+    params = [
+        (vmat0, 2, 0.5),
+        (vmat0, 4, 0.5),
+        (vmat1, 2, op_norm(2)),
+        (vmat1, 4, op_norm(4)),
+        (vmat2, 2, 6 * op_norm(2) ** 2),
+        (vmat2, 4, 6 * op_norm(4) ** 2),
+        (vmat3, 2, np.sqrt(18 * op_norm(2) ** 3)),
+        (vmat3, 4, np.sqrt(18 * op_norm(4) ** 3)),
+        (vmat4, 2, 6 * op_norm(2) ** 2 + np.sqrt(18 * op_norm(2) ** 3)),
+        (vmat4, 4, 6 * op_norm(4) ** 2 + np.sqrt(18 * op_norm(4) ** 3)),
+    ]
+
+    @pytest.mark.parametrize("vmatrix, gridpoints, expected", params)
+    def test_norm(self, vmatrix: VibronicMatrix, gridpoints: int, expected: float):
+        """Test that the norm is correct"""
+
+        assert np.isclose(vmatrix.norm(gridpoints), expected)
+
+    params = [
+        (vmat0, 2),
+        (vmat0, 4),
+        (vmat1, 2),
+        (vmat1, 4),
+        (vmat2, 2),
+        (vmat2, 4),
+        (vmat3, 2),
+        (vmat3, 4),
+        (vmat4, 2),
+        (vmat4, 4),
+    ]
+
+    @pytest.mark.parametrize("vmatrix, gridpoints", params)
+    def test_norm_against_numpy(self, vmatrix: VibronicMatrix, gridpoints: int):
+        """Test that .norm is an upper bound on the true norm"""
+        upper_bound = vmatrix.norm(gridpoints)
+        norm = np.max(np.linalg.eigvals(vmatrix.matrix(gridpoints)))
+        assert norm <= upper_bound