perf(numba): drop the redundant passes in curvature-matrix F assembly (#505)

Three redundant passes over F in the numba sparse imaging inversion, none of
which changed a single value:

1. A global `curvature_matrix_mirrored_from` ran over the whole (P, P) matrix
   after assembly. Every block was already symmetric or had a known transpose:
   the mapper x mapper blocks are folded and mirrored inside
   `curvature_matrix_via_sparse_operator_from`, and the linear-func x
   linear-func blocks already wrote both triangles. The mirror is removed and
   the two off-diagonal block writers (mapper x mapper, mapper x linear-func)
   now place their transpose alongside the block, so F leaves assembly
   symmetric.
2. `_curvature_matrix_mapper_diag` wrapped the three `mapper.unique_mappings`
   arrays in `np.array(...)`, copying them on every evaluation. They are
   already contiguous ndarrays of the dtype the numba kernel wants (see
   `UniqueMappings.__init__`, which casts on construction) and every other
   caller in this module passes them through untouched.
3. `_curvature_matrix_mapper_func_blocks_from` formed
   `operated_mapping_matrix / noise_map ** 2` inside the mapper loop and then
   copied the result again with `np.array`. The weights do not depend on the
   mapper, so they are formed once per linear func ahead of the loop and passed
   straight to the kernel.

`curvature_matrix_mirrored_from` in `inversion_imaging_numba_util` had no
remaining caller (in this repo or downstream) and is deleted rather than left
as dead numba code.

Also corrects the `curvature_matrix` docstring: it claimed the property is
"not a cached property" and is overwritten in memory by the regularization
add, but it is decorated `@cached_property` and `curvature_reg_matrix` adds
out-of-place via `np.add`.

Bit-identical output verified by capturing `inversion.curvature_matrix` before
and after on the autolens_profiling breakdown fiducial: `np.array_equal` True
for both euclid and hst, with the figures of merit matching the pins exactly
(hst 27661.910133664103, euclid 6213.3068738858765).

`pytest test_autoarray/inversion`: 392 passed.

Step 1 of #505.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Fr6iJ5T1RDARWfxttCuGkK

Author: Jammy2211

autoarray/inversion/inversion/imaging_numba/inversion_imaging_numba_util.py

@@ -492,26 +492,6 @@ def curvature_matrix_with_added_to_diag_from(
     return curvature_matrix
 
 
-@numba_util.jit()
-def curvature_matrix_mirrored_from(
-    curvature_matrix: np.ndarray,
-) -> np.ndarray:
-    curvature_matrix_mirrored = np.zeros(
-        (curvature_matrix.shape[0], curvature_matrix.shape[1])
-    )
-
-    for i in range(curvature_matrix.shape[0]):
-        for j in range(curvature_matrix.shape[1]):
-            if curvature_matrix[i, j] != 0:
-                curvature_matrix_mirrored[i, j] = curvature_matrix[i, j]
-                curvature_matrix_mirrored[j, i] = curvature_matrix[i, j]
-            if curvature_matrix[j, i] != 0:
-                curvature_matrix_mirrored[i, j] = curvature_matrix[j, i]
-                curvature_matrix_mirrored[j, i] = curvature_matrix[j, i]
-
-    return curvature_matrix_mirrored
-
-
 @numba_util.jit()
 def curvature_matrix_via_sparse_operator_from(
     psf_precision_operator: np.ndarray,

autoarray/inversion/inversion/imaging_numba/sparse.py

@@ -332,10 +332,14 @@ def curvature_matrix(self) -> np.ndarray:
         for simultaneously. In the w-tilde formalism this requires us to consider the mappings between data and every
         linear object, meaning that the linear alegbra has both on and off diagonal terms.
 
-        The `curvature_matrix` computed here is overwritten in memory when the regularization matrix is added to it,
-        because for large matrices this avoids overhead. For this reason, `curvature_matrix` is not a cached property
-        to ensure if we access it after computing the `curvature_reg_matrix` it is correctly recalculated in a new
-        array of memory.
+        Every block of F is written into the matrix already symmetrized: the mapper x mapper blocks are
+        folded and mirrored inside `curvature_matrix_via_sparse_operator_from`, and the off-diagonal
+        blocks (mapper x mapper, mapper x linear-func, linear-func x linear-func) are each placed
+        together with their transpose. A global symmetrizing pass over the assembled matrix would
+        therefore be a no-op, and is not run.
+
+        `curvature_matrix` is a cached property, and `curvature_reg_matrix` adds the regularization
+        matrix to it out-of-place, so the cached F is never overwritten by that addition.
         """
         if self.has(cls=AbstractLinearObjFuncList):
             curvature_matrix = self._curvature_matrix_func_list_and_mapper
@@ -344,10 +348,6 @@ def curvature_matrix(self) -> np.ndarray:
         else:
             curvature_matrix = self._curvature_matrix_multi_mapper
 
-        curvature_matrix = inversion_imaging_numba_util.curvature_matrix_mirrored_from(
-            curvature_matrix=curvature_matrix,
-        )
-
         if len(self.no_regularization_index_list) > 0:
             curvature_matrix = (
                 inversion_imaging_numba_util.curvature_matrix_with_added_to_diag_from(
@@ -386,11 +386,9 @@ def _curvature_matrix_mapper_diag(self) -> Optional[np.ndarray]:
                 psf_precision_operator=self.sparse_operator.psf_precision_operator_sparse,
                 psf_precision_indexes=self.sparse_operator.indexes,
                 psf_precision_lengths=self.sparse_operator.lengths,
-                data_to_pix_unique=np.array(
-                    mapper_i.unique_mappings.data_to_pix_unique
-                ),
-                data_weights=np.array(mapper_i.unique_mappings.data_weights),
-                pix_lengths=np.array(mapper_i.unique_mappings.pix_lengths),
+                data_to_pix_unique=mapper_i.unique_mappings.data_to_pix_unique,
+                data_weights=mapper_i.unique_mappings.data_weights,
+                pix_lengths=mapper_i.unique_mappings.pix_lengths,
                 pix_pixels=mapper_i.params,
             )
 
@@ -493,6 +491,11 @@ def _curvature_matrix_multi_mapper(self) -> np.ndarray:
                     mapper_param_range_j[0] : mapper_param_range_j[1],
                 ] = off_diag
 
+                curvature_matrix[
+                    mapper_param_range_j[0] : mapper_param_range_j[1],
+                    mapper_param_range_i[0] : mapper_param_range_i[1],
+                ] = off_diag.T
+
         return curvature_matrix
 
     @property
@@ -540,8 +543,9 @@ def _curvature_matrix_mapper_func_blocks_from(
         noise-weighted curvature vector of a linear function (the dense sliding-window correlation in
         `curvature_matrix_off_diags_via_mapper_and_linear_func_curvature_vector_from`).
 
-        Only the `[mapper, linear_func]` blocks are written; their transposes are filled in by the
-        global mirror applied in `curvature_matrix`.
+        Each `[mapper, linear_func]` block is written together with its transpose into the
+        `[linear_func, mapper]` block, so F leaves this helper symmetric and no global mirroring
+        pass is required.
 
         Parameters
         ----------
@@ -557,24 +561,29 @@ def _curvature_matrix_mapper_func_blocks_from(
             cls=AbstractLinearObjFuncList
         )
 
+        # The noise-weighted curvature weights of a linear func do not depend on the mapper, so
+        # they are formed once per linear func rather than once per (mapper, linear func) pair.
+        curvature_weights_list = [
+            np.array(
+                self.linear_func_operated_mapping_matrix_dict[linear_func]
+                / self.noise_map[:, None] ** 2
+            )
+            for linear_func in linear_func_list
+        ]
+
         for i in range(len(mapper_list)):
             mapper = mapper_list[i]
             mapper_param_range = mapper_param_range_list[i]
 
             for func_index, linear_func in enumerate(linear_func_list):
                 linear_func_param_range = linear_func_param_range_list[func_index]
 
-                data_linear_func_matrix = (
-                    self.linear_func_operated_mapping_matrix_dict[linear_func]
-                    / self.noise_map[:, None] ** 2
-                )
-
                 off_diag = inversion_imaging_numba_util.curvature_matrix_off_diags_via_mapper_and_linear_func_curvature_vector_from(
                     data_to_pix_unique=mapper.unique_mappings.data_to_pix_unique,
                     data_weights=mapper.unique_mappings.data_weights,
                     pix_lengths=mapper.unique_mappings.pix_lengths,
                     pix_pixels=mapper.params,
-                    curvature_weights=np.array(data_linear_func_matrix),
+                    curvature_weights=curvature_weights_list[func_index],
                     mask=self.mask.array,
                     psf_kernel=self.psf.kernel.native.array,
                 )
@@ -584,6 +593,11 @@ def _curvature_matrix_mapper_func_blocks_from(
                     linear_func_param_range[0] : linear_func_param_range[1],
                 ] = off_diag
 
+                curvature_matrix[
+                    linear_func_param_range[0] : linear_func_param_range[1],
+                    mapper_param_range[0] : mapper_param_range[1],
+                ] = off_diag.T
+
         return curvature_matrix
 
     def _curvature_matrix_func_func_blocks_from(