Speed improvements to inversion for large datasets

[billetd]

These changes make improvements to the speed of inversions when there is a lot of data.

I was running into slow inversion issues(?) in my use of Lompe for widebeam SuperDARN data. The inversion process (model.run_inversion()) was by far the slowest part of the process, mainly because there was so much data in the inversion, making the SVD scipy solver alone take almost a full second per inversion.

This PR makes two changes:

• Replace the SVD part with a Cholesky solver, with a fallback to SVD. Cholesky should produce mathematically identical results as long as the array is non-zero. SVD still handles edge cases.
• Refactors how GTG_i and GTD_i so it does it via broadcasting rather than numpy math. Avoids constantly remaking the matrix.

I have tested this on my FBI library for high-res SuperDARN measurements and noticed about a ~50% improvement in inversion speeds (~0.5 vs ~1s) from the previous formulation with bitwise-identical results. I have also tested in a more general Lompe use-case using an adapted version of two_ways_to_calculate_FAC.py:

import lompe
import numpy as np
import pandas as pd
from datetime import datetime
from lompe import lompeplot
from lompe.utils.conductance import hardy_EUV
import apexpy
import matplotlib.pyplot as plt
import time

event = '2012-04-05'
apex = apexpy.Apex(int(event[0:4]), refh = 110)

supermagfn = '/Users/danielbillett/PycharmProjects/lompe/examples/sample_dataset/20120405_supermag.h5'
superdarnfn = '/Users/danielbillett/PycharmProjects/lompe/examples/sample_dataset/20120405_superdarn_grdmap.h5'
iridiumfn = '/Users/danielbillett/PycharmProjects/lompe/examples/sample_dataset/20120405_iridium.h5'

# set up grid
position = (-90, 65) # lon, lat
orientation = (-1, 2) # east, north
L, W, Lres, Wres = 4200e3, 7000e3, 100.e3, 100.e3 # dimensions and resolution of grid
grid = lompe.cs.CSgrid(lompe.cs.CSprojection(position, orientation), L, W, Lres, Wres, R = 6481.2e3)

# load ampere, supermag, and superdarn data from 2012-05-05
ampere    = pd.read_hdf(iridiumfn)
supermag  = pd.read_hdf(supermagfn)
superdarn = pd.read_hdf(superdarnfn)

t0, t1 = datetime(2012, 4, 5, 12, 0), datetime(2012, 4, 5, 12, 5)

# crop AMPERE data and make lompe Data object:
amp = ampere[(ampere.time >= t0) & (ampere.time <= t1)]
B = np.vstack((amp.B_e.values, amp.B_n.values, amp.B_r.values))
coords = np.vstack((amp.lon.values, amp.lat.values, amp.r.values))
amp_data = lompe.Data(B * 1e-9, coords, datatype = 'space_mag_fac')

# crop supermag and make lompe Data object
sm = supermag[t0:t1]
B = np.vstack((sm.Be.values, sm.Bn.values, sm.Bu.values))
coords = np.vstack((sm.lon.values, sm.lat.values))
sm_data = lompe.Data(B * 1e-9, coords, datatype = 'ground_mag')

# crop superdarn and make lompe Data object
sd = superdarn.loc[(superdarn.index >= t0) & (superdarn.index <= t1)]
vlos = sd['vlos'].values
coords = np.vstack((sd['glon'].values, sd['glat'].values))
los  = np.vstack((sd['le'].values, sd['ln'].values))
sd_data = lompe.Data(vlos, coordinates = coords, LOS = los, datatype = 'convection')


Kp = 4 # for Hardy conductance model
SH = lambda lon = grid.lon, lat = grid.lat: hardy_EUV(lon, lat, 5, t0, 'hall'    )
SP = lambda lon = grid.lon, lat = grid.lat: hardy_EUV(lon, lat, 5, t0, 'pedersen')
model = lompe.Emodel(grid, Hall_Pedersen_conductance = (SH, SP))

start_time = time.perf_counter()
model.add_data(amp_data, sm_data, sd_data)
end_time = time.perf_counter()
execution_time = end_time - start_time
print(f"Execution time: {execution_time} seconds")
model.run_inversion(l1 = 1, l2 = 10)


lompeplot(model, show_data_location=True)
plt.show()

This also gave identical results, albeit the speed improvement was more like ~10% due to the lower amount of data involved.

I do think this needs to be tested well if it's ever merged. I'm pretty confident it's a free speedup, but there could potentially be edge cases where a Cholesky solver gives different results to SVD (e.g., if GG is wholly positive, but has very small values in it. I haven't found a case yet where the results aren't identical.

[BingMM]

A couple of comments. Strictly speaking, the requirement is that the system matrix is symmetric positive definite (SPD), not merely non-zero. Cholesky factorization can still succeed for an extremely ill-conditioned matrix, so success alone does not guarantee a well-conditioned solve. One could perform additional diagnostics on the factorization, but that is probably overkill here and would work against the goal of improving performance.

An additional speed-up could be achieved by bypassing the computation of the posterior covariance. In most use cases, only the MAP/model estimate is required. As long as the Cholesky factorization is retained, the posterior covariance can be computed later if and when it is needed.

This is essentially the approach used in nandir