code: def _compress_rank(white_vec, prec_sqrt, assume_full_rank=False):
      """
      Compress a wide representation ``(white_vec, prec_sqrt)`` while preserving
      the quadratic function ``||x @ prec_sqrt - white_vec||^2 + const``.
      """
      dim, rank = prec_sqrt.shape[-2:]
      assert rank >= dim
      old_norm2 = _norm2(white_vec)
  
      # Let P = prec_sqrt and w = white_vec define the original Gaussian
      #
      #   G(x;w,P) = -1/2 || x P - w ||^2
      #            = -1/2 x P P' x' + x P w' -1/2 w w'
      #
      # We seek a compressed Gaussian G(x;wc,Pc) and constant C such that
      #
      #   G(x;w,P) = G(x;wc,Pc) + C
      #            = -1/2 x Pc Pc' x' + x Pc wc' -1/2 wc wc' + C
      if assume_full_rank:
          # Cholesky factorizing Pc = chol(P P'), we match remaining coefficients
          #
          #    Pc wc' = P w'  ==>  wc' = Pc \ P w'
          info_vec_ = prec_sqrt @ white_vec[..., None]
          precision = prec_sqrt @ ops.transpose(prec_sqrt, -1, -2)
          prec_sqrt = ops.cholesky(precision)
          white_vec = ops.triangular_solve(info_vec_, prec_sqrt)[..., 0]
      else:
          # Computing a reduced QR representation of P' of shape (rank,dim)
          #
          #   P' = Q [ R ]     P = [ R'  0 ] Q'
          #          [ 0 ]
          #
          # where Q is orthogonal and R is upper triangular of shape (dim,dim).
          # Then splitting along the new dimension,
          #
          #   G(x;w,P) = -1/2 || x [R' 0] Q' - w ||^2
          #            = -1/2 || x [R' 0] - w Q ||^2
          #            = -1/2 || x R' - (w Q)[...,:dim] ||^2
          #              -1/2 || (w Q)[...,dim:] ||^2
          #            =: G(x;wc,Pc) + C
          # where
          #
          #   wc = (w Q)[...,:dim]
          #   Pc = R'
          Q, R = ops.qr(ops.transpose(prec_sqrt, -1, -2))
          assert Q.shape[-2:] == (
              rank,
              dim,
          )  # note only part of Q is returned
          assert R.shape[-2:] == (dim, dim)
          prec_sqrt = ops.transpose(R, -1, -2)
          white_vec = _vm(white_vec, Q)
      # Finally the shifting constant is
      #
      #    C = 1/2 (wc wc' - w w')
      new_norm2 = _norm2(white_vec)
      shift = 0.5 * (new_norm2 - old_norm2)
      return white_vec, prec_sqrt, shift
  
imports: import math
  from collections import OrderedDict, defaultdict
  from contextlib import contextmanager
  from functools import reduce
  import funsor.ops
  from funsor.affine import affine_inputs, extract_affine, is_affine
  from funsor.delta import Delta
  from funsor.domains import Real, Reals
  from funsor.interpretations import compress_gaussians
  from funsor.ops import AddOp, SubOp
  from funsor.tensor import Tensor, align_tensor, align_tensors
  from funsor.terms import Binary, Funsor, FunsorMeta, Number, Slice, Subs, Variable, eager, reflect
  from funsor.util import broadcast_shape, get_tracing_state, lazy_property
label_desc: Use image node_id 9 for calling a global function with example usage: _compress_rank(white_vec, prec_sqrt, assume_full_rank) and returns: white_vec,  prec_sqrt,  shift