integrated PDMats.jl for speedup

[riteshbhirud]

PDMats.jl Integration for Efficient Covariance Matrix Operations

Fixes #53

Summary

This pull request integrates PDMats.jl into Gabs.jl to optimize operations involving covariance matrices in Gaussian states. By leveraging specialized positive definite matrix types (e.g. ScalMat, PDiagMat, and PDMat), this integration achieves significant performance improvements without sacrificing numerical accuracy or API consistency.

Performance Improvements

Benchmarking demonstrates dramatic speedups for key operations as the system size increases:
• Inverse quadratic forms: up to 490× faster for 20-mode systems
• Whitening operations: up to 233× faster for 20-mode systems
• State creation: up to 82× faster for 20-mode systems
• Purity calculation: up to 9.4× faster for 20-mode systems
• Partial trace operations: approximately 2× faster than the standard implementation

Implementation Details

New PDGaussianState Type

A new type, PDGaussianState, extends the standard GaussianState by storing the covariance matrix as an optimized PDMats type (AbstractPDMat). This allows efficient linear algebra routines to be used for operations such as:
• Purity calculation
• Whitening and unwhitening
• Quadratic and inverse quadratic forms

Optimized Constructors

New constructors have been added to build PDMats‑optimized states:
• vacuumstate_pd
• thermalstate_pd
• coherentstate_pd

These constructors ensure that the covariance matrices are constructed directly as PDMats types (e.g. ScalMat or PDiagMat when appropriate) for maximum efficiency.

Core Operations

All major operations have been updated to use the PDGaussianState covariance:
• Purity calculation, whitening/unwhitening, and quadratic forms: These now call the corresponding PDMats routines.
• Tensor product: A specialized tensor routine builds the block–diagonal covariance directly from the PD matrices using a helper function that dispatches on the covariance type.
• Partial trace: New optimized implementations extract the appropriate 2×2 block for a given mode (or modes) directly from the underlying PD covariance—without converting the entire matrix to dense form—thus providing substantial performance improvements.

Usage Examples

The API remains consistent with the standard GaussianState, so existing code will work seamlessly with the new PD-optimized functions. For example:`using Gabs

Create a standard 10-mode state

basis = QuadPairBasis(10)
std_state = vacuumstate(basis)

Create a PDMats-optimized 10-mode state

pd_state = vacuumstate_pd(basis)

Perform operations using the same API:

x = randn(2 * basis.nmodes)
whitened_x = whiten(pd_state, x)
quad_value = invquad(pd_state, x)

Tensor product and partial trace:

tensor_state = tensor(pd_state, pd_state)
reduced_state = ptrace(tensor_state, 1)

println("Whitening result: ", whitened_x)
println("Inverse quadratic form: ", quad_value)
println("Partial trace mean: ", reduced_state.mean)
println("Partial trace covariance: ", reduced_state.covar)`

Numerical Validation

A comprehensive test suite confirms that the new PDGaussianState operations (including the optimized partial trace) produce results that are numerically equivalent to the original implementation. Benchmarks and unit tests verify that the improvements in speed come with exact accuracy, ensuring full functionality.

[apkille]

I haven't looked at everything in great detail yet, but here are a few things:

• Why do you have to define a whole new composite type for PDGaussianState? It has the same fields as GaussianState right? So don't we just need to define functions and constants similar to SparseOperator in QuantumOptics.jl?
• Rather than defining helper functions like name_pd, lets just the dispatch method for calling PDGaussianStateType types as documented here.
• Let's add the benchmarks that you put to the src/benchmarks folder rather than adding it as a dependency and into the tests.

Let me know if you have any questions!

[apkille]

Marking this as a draft to make my review log cleaner. When it's ready please ping me!