diff --git a/.gitignore b/.gitignore index e228a4bb..4c8fb780 100644 --- a/.gitignore +++ b/.gitignore @@ -17,3 +17,4 @@ Cargo.lock # coverage *.xml +ndarray-linalg/build.rs diff --git a/ndarray-linalg/Cargo.toml b/ndarray-linalg/Cargo.toml index e06fa675..91edb473 100644 --- a/ndarray-linalg/Cargo.toml +++ b/ndarray-linalg/Cargo.toml @@ -35,6 +35,7 @@ num-complex = "0.4.0" num-traits = "0.2.14" rand = "0.8.3" thiserror = "1.0.24" +statrs = "0.16.0" [dependencies.ndarray] version = "0.15.2" diff --git a/ndarray-linalg/src/expm.rs b/ndarray-linalg/src/expm.rs new file mode 100644 index 00000000..47267e6e --- /dev/null +++ b/ndarray-linalg/src/expm.rs @@ -0,0 +1,297 @@ +use crate::{Inverse, OperationNorm}; +use cauchy::Scalar; +use lax::Lapack; +use ndarray::prelude::*; +extern crate statrs; +use crate::error::{LinalgError, Result}; +use statrs::function::factorial::{binomial, factorial}; + +// These constants are hard-coded from Al-Mohy & Higham +const THETA_3: f64 = 1.495_585_217_958_292e-2; +const THETA_5: f64 = 2.539_398_330_063_23e-1; + +const THETA_7: f64 = 9.504_178_996_162_932e-1; +const THETA_9: f64 = 2.097_847_961_257_068e0; +const THETA_13: f64 = 4.25; // Alg 5.1 + +// The Pade Coefficients for the numerator of the diagonal approximation to Exp[x]. Computed via Mathematica/WolframAlpha. +// Note that the denominator has the same coefficients but odd powers have an opposite sign. +// Coefficients are also stored via the power of x, so for example the numerator would look like +// x^0 * PADE_COEFF_M[0] + x^1 * PADE_COEFF_M[1] + ... + x^m * PADE_COEFF_M[M] +const PADE_COEFFS_3: [f64; 4] = [1., 0.5, 0.1, 1. / 120.]; + +const PADE_COEFFS_5: [f64; 6] = [1., 0.5, 1. / 9., 1. / 72., 1. / 1_008., 1. / 30_240.]; + +const PADE_COEFFS_7: [f64; 8] = [ + 1., + 0.5, + 3. / 26., + 5. / 312., + 5. / 3_432., + 1. / 11_440., + 1. / 308_880., + 1. / 17_297_280., +]; + +const PADE_COEFFS_9: [f64; 10] = [ + 1., + 0.5, + 2. / 17., + 7. / 408., + 7. / 4_080., + 1. / 8_160., + 1. / 159_120., + 1. / 4_455_360., + 1. / 196_035_840., + 1. / 17_643_225_600., +]; + +const PADE_COEFFS_13: [f64; 14] = [ + 1., + 0.5, + 3. / 25., + 11. / 600., + 11. / 5_520., + 3. / 18_400., + 1. / 96_600., + 1. / 1_932_000., + 1. / 48_944_000., + 1. / 1_585_785_600., + 1. / 67_395_888_000., + 1. / 3_953_892_096_000., + 1. / 355_850_288_640_000., + 1. / 64_764_752_532_480_000., +]; + +fn pade_approximation_3 + Lapack>( + a_1: &Array2, + a_2: &Array2, +) -> Result> { + let mut evens: Array2 = Array2::::eye(a_1.nrows()); + // evens.mapv_inplace(|x| x * S::from_real(PADE_COEFFS_3[0])); + evens.scaled_add(S::from_real(PADE_COEFFS_3[2]), a_2); + + let mut odds: Array2 = Array2::::eye(a_1.nrows()); + odds.mapv_inplace(|x| x * S::from_real(PADE_COEFFS_3[1])); + odds.scaled_add(S::from_real(PADE_COEFFS_3[3]), a_2); + odds = odds.dot(a_1); + + odds.mapv_inplace(|x| -x); + let inverted = (&odds + &evens).inv()?; + odds.mapv_inplace(|x| -x); + Ok(inverted.dot(&(odds + evens))) +} + +fn pade_approximation_5 + Lapack>( + a_1: &Array2, + a_2: &Array2, + a_4: &Array2, +) -> Result> { + let mut evens: Array2 = Array2::::eye(a_1.nrows()); + // evens.mapv_inplace(|x| S::from_real(PADE_COEFFS_5[0]) * x); + evens.scaled_add(S::from_real(PADE_COEFFS_5[2]), a_2); + evens.scaled_add(S::from_real(PADE_COEFFS_5[4]), a_4); + + let mut odds: Array2 = Array::eye(a_1.nrows()); + odds.mapv_inplace(|x| S::from_real(PADE_COEFFS_5[1]) * x); + odds.scaled_add(S::from_real(PADE_COEFFS_5[3]), a_2); + odds.scaled_add(S::from_real(PADE_COEFFS_5[5]), a_4); + odds = odds.dot(a_1); + + odds.mapv_inplace(|x| -x); + let inverted = (&odds + &evens).inv()?; + odds.mapv_inplace(|x| -x); + Ok(inverted.dot(&(odds + evens))) +} + +fn pade_approximation_7 + Lapack>( + a_1: &Array2, + a_2: &Array2, + a_4: &Array2, + a_6: &Array2, +) -> Result> { + let mut evens: Array2 = Array::eye(a_1.nrows()); + // evens.mapv_inplace(|x| S::from_real(PADE_COEFFS_7[0]) * x); + evens.scaled_add(S::from_real(PADE_COEFFS_7[2]), a_2); + evens.scaled_add(S::from_real(PADE_COEFFS_7[4]), a_4); + evens.scaled_add(S::from_real(PADE_COEFFS_7[6]), a_6); + + let mut odds: Array2 = Array::eye(a_1.nrows()); + odds.mapv_inplace(|x| S::from_real(PADE_COEFFS_7[1]) * x); + odds.scaled_add(S::from_real(PADE_COEFFS_7[3]), a_2); + odds.scaled_add(S::from_real(PADE_COEFFS_7[5]), a_4); + odds.scaled_add(S::from_real(PADE_COEFFS_7[7]), a_6); + odds = odds.dot(a_1); + + odds.mapv_inplace(|x| -x); + let inverted = (&odds + &evens).inv()?; + odds.mapv_inplace(|x| -x); + Ok(inverted.dot(&(odds + evens))) +} + +fn pade_approximation_9 + Lapack>( + a_1: &Array2, + a_2: &Array2, + a_4: &Array2, + a_6: &Array2, + a_8: &Array2, +) -> Result> { + let mut evens: Array2 = Array::eye(a_1.nrows()); + // evens.mapv_inplace(|x| S::from_real(PADE_COEFFS_9[0]) * x); + evens.scaled_add(S::from_real(PADE_COEFFS_9[2]), a_2); + evens.scaled_add(S::from_real(PADE_COEFFS_9[4]), a_4); + evens.scaled_add(S::from_real(PADE_COEFFS_9[6]), a_6); + evens.scaled_add(S::from_real(PADE_COEFFS_9[8]), a_8); + + let mut odds: Array2 = Array::eye(a_1.nrows()); + odds.mapv_inplace(|x| S::from_real(PADE_COEFFS_9[1]) * x); + odds.scaled_add(S::from_real(PADE_COEFFS_9[3]), a_2); + odds.scaled_add(S::from_real(PADE_COEFFS_9[5]), a_4); + odds.scaled_add(S::from_real(PADE_COEFFS_9[7]), a_6); + odds.scaled_add(S::from_real(PADE_COEFFS_9[9]), a_8); + odds = odds.dot(a_1); + + odds.mapv_inplace(|x| -x); + let inverted: Array2 = (&odds + &evens).inv()?; + odds.mapv_inplace(|x| -x); + Ok(inverted.dot(&(odds + evens))) +} + +// Note: all input matrices should be scaled in the main expm +// function. +fn pade_approximation_13 + Lapack>( + a_1: &Array2, + a_2: &Array2, + a_4: &Array2, + a_6: &Array2, +) -> Result> { + let mut evens_1: Array2 = Array::eye(a_1.nrows()); + evens_1.mapv_inplace(|x| S::from_real(PADE_COEFFS_13[0]) * x); + evens_1.scaled_add(S::from_real(PADE_COEFFS_13[2]), a_2); + evens_1.scaled_add(S::from_real(PADE_COEFFS_13[4]), a_4); + evens_1.scaled_add(S::from_real(PADE_COEFFS_13[6]), a_6); + + let mut evens_2 = a_2.clone(); + evens_2.mapv_inplace(|x| S::from_real(PADE_COEFFS_13[8]) * x); + evens_2.scaled_add(S::from_real(PADE_COEFFS_13[10]), a_4); + evens_2.scaled_add(S::from_real(PADE_COEFFS_13[12]), a_6); + let evens = evens_2.dot(a_6) + &evens_1; + + let mut odds_1: Array2 = Array::eye(a_1.nrows()); + odds_1.mapv_inplace(|x| S::from_real(PADE_COEFFS_13[1]) * x); + odds_1.scaled_add(S::from_real(PADE_COEFFS_13[3]), a_2); + odds_1.scaled_add(S::from_real(PADE_COEFFS_13[5]), a_4); + odds_1.scaled_add(S::from_real(PADE_COEFFS_13[7]), a_6); + + let mut odds_2 = a_2.clone(); + odds_2.mapv_inplace(|x| S::from_real(PADE_COEFFS_13[9]) * x); + odds_2.scaled_add(S::from_real(PADE_COEFFS_13[11]), a_4); + odds_2.scaled_add(S::from_real(PADE_COEFFS_13[13]), a_6); + odds_2 = odds_2.dot(a_6); + + let mut odds = (&odds_1 + &odds_2).dot(a_1); + odds.mapv_inplace(|x| -x); + let inverted: Array2 = (&odds + &evens).inv()?; + odds.mapv_inplace(|x| -x); + Ok(inverted.dot(&(odds + evens))) +} + +fn power_abs_norm(input_matrix: &Array2, p: usize) -> f64 +where + S: Scalar, +{ + let mut v = Array1::::ones((input_matrix.ncols()).f()); + let abs_matrix = input_matrix.t().map(|x| x.abs()); + for _ in 0..p { + v.assign(&abs_matrix.dot(&v)); + } + // return max col sum + v.into_iter() + .reduce(|x, y| if x > y { x } else { y }) + .unwrap() +} + +/// helper function used in Al-Mohy & Higham. Name is unchanged for both literature reference. +fn ell>(a_matrix: &Array2, m: u64) -> Result { + if a_matrix.is_square() == false { + return Err(LinalgError::NotSquare { + rows: a_matrix.nrows() as i32, + cols: a_matrix.ncols() as i32, + }); + } + let p = 2 * m + 1; + let a_one_norm = a_matrix.map(|x| x.abs()).opnorm_one()?; + + let powered_abs_norm = power_abs_norm(a_matrix, p as usize); + let alpha = powered_abs_norm * pade_error_coefficient(m) / a_one_norm; + let u = f64::EPSILON / 2.; + let log2_alpha_div_u = f64::log2(alpha / u); + let val = f64::ceil(log2_alpha_div_u / ((2 * m) as f64)) as i32; + Ok(i32::max(val, 0)) +} + +/// Calculates the leading term of the error series for the [m/m] Pade approximation to exp(x). +fn pade_error_coefficient(m: u64) -> f64 { + 1.0 / (binomial(2 * m, m) * factorial(2 * m + 1)) +} + +/// ## Matrix Exponentiation +/// Computes matrix exponential based on the scaling-and-squaring algorithm by Al-Mohy and Higham [[1]]. +/// Currently restricted to matrices with entries that are either f64 or Complex64. 64 bit precision is required +/// due to error calculations in [[1]] Utilizes Lapack calls so entries must satisfy LAPACK trait bounds. +/// +/// ### Caveats +/// Currently confirmed accurate to f64 precision up to 1024x1024 sparse matrices. Dense matrices +/// have been observed with a worse average entry error, up to 100x100 matrices should be close +/// enough to f64 precision for vast majority of numeric purposes. +/// +/// ### References +/// [[1]] A New Scaling and Squaring Algorithm for the Matrix Exponential. Al-Mohy, Awad H. and Higham, Nicholas J. 2009. Siam J. Matrix Anal. Appl. Vol. 31, No. 3, pp. 970-989. +pub fn expm + Lapack>(a_matrix: &Array2) -> Result> { + let mut a_2 = a_matrix.dot(a_matrix); + let mut a_4 = a_2.dot(&a_2); + let mut a_6 = a_2.dot(&a_4); + let d4 = a_4.opnorm_one()?.powf(1. / 4.); + let d6 = a_6.opnorm_one()?.powf(1. / 6.); + // Note d6 should be an estimate and d4 an estimate + let eta_1 = f64::max(d4, d6); + if eta_1 < THETA_3 && ell(a_matrix, 3)? == 0 { + return pade_approximation_3(a_matrix, &a_2); + } + // d4 should be exact here, d6 an estimate + let eta_2 = f64::max(d4, d6); + if eta_2 < THETA_5 && ell(a_matrix, 5)? == 0 { + return pade_approximation_5(a_matrix, &a_2, &a_4); + } + let a_8 = a_4.dot(&a_4); + let d8 = a_8.opnorm_one().unwrap().powf(1. / 8.); + let eta_3 = f64::max(d6, d8); + if eta_3 < THETA_7 && ell(a_matrix, 7)? == 0 { + return pade_approximation_7(a_matrix, &a_2, &a_4, &a_6); + } + if eta_3 < THETA_9 && ell(a_matrix, 9)? == 0 { + return pade_approximation_9(a_matrix, &a_2, &a_4, &a_6, &a_8); + } + let a_10 = a_2.dot(&a_8); + let eta_4 = f64::max(d8, a_10.opnorm_one()?); + let eta_5 = f64::min(eta_3, eta_4); + + let mut s = f64::max(0., (eta_5 / THETA_13).log2().ceil()) as i32; + let mut a_scaled = a_matrix.clone(); + let mut scaler = S::from_real(2.).powi(-s); + a_scaled.mapv_inplace(|x| x * scaler); + s += ell(&a_scaled, 13)?; + + a_scaled.assign(a_matrix); + scaler = S::from_real(2.).powi(-s); + a_scaled.mapv_inplace(|x| x * scaler); + a_2.mapv_inplace(|x| x * scaler.powi(2)); + a_4.mapv_inplace(|x| x * scaler.powi(4)); + a_6.mapv_inplace(|x| x * scaler.powi(6)); + + let mut output = pade_approximation_13(&a_scaled, &a_2, &a_4, &a_6)?; + for _ in 0..s { + output = output.dot(&output); + } + Ok(output) +} diff --git a/ndarray-linalg/src/lib.rs b/ndarray-linalg/src/lib.rs index 784e1dff..bd6d609e 100644 --- a/ndarray-linalg/src/lib.rs +++ b/ndarray-linalg/src/lib.rs @@ -56,6 +56,7 @@ pub mod diagonal; pub mod eig; pub mod eigh; pub mod error; +pub mod expm; pub mod generate; pub mod inner; pub mod krylov; @@ -81,6 +82,7 @@ pub use crate::convert::*; pub use crate::diagonal::*; pub use crate::eig::*; pub use crate::eigh::*; +pub use crate::expm::expm; pub use crate::generate::*; pub use crate::inner::*; pub use crate::layout::*; diff --git a/ndarray-linalg/tests/expm.rs b/ndarray-linalg/tests/expm.rs new file mode 100644 index 00000000..37ef30c1 --- /dev/null +++ b/ndarray-linalg/tests/expm.rs @@ -0,0 +1,109 @@ +use ndarray::linalg::kron; +use ndarray::*; +use ndarray_linalg::expm::expm; +use ndarray_linalg::{Eig, OperationNorm}; +use num_complex::{Complex64 as c64, ComplexFloat}; +use rand::*; + +/// Test matrix exponentiation expm by using an exact formula for exponentials of Pauli matrices. +/// These matrices are 1-sparse which allows for a different testing regime than the dense matrix +/// test. +#[test] +fn test_random_sparse_matrix() { + let mut rng = rand::thread_rng(); + let num_qubits: u32 = 10; + let dim = 2_usize.pow(num_qubits); + let mut matrix = Array2::::eye(2); + let _i = c64::new(0., 1.); + let zero = c64::new(0., 0.); + let pauli_x = array![[zero, c64::new(1., 0.)], [c64::new(1., 0.), zero]]; + let pauli_y = array![[zero, c64::new(0., -1.)], [c64::new(0., 1.), zero]]; + let pauli_z = array![[c64::new(1., 0.), zero], [zero, c64::new(-1., 0.)]]; + for n in 0..num_qubits { + let pauli_matrix = match rng.gen_range::(0..=3) { + 0 => Array2::::eye(2), + 1 => pauli_x.clone(), + 2 => pauli_y.clone(), + 3 => pauli_z.clone(), + _ => unreachable!(), + }; + if n == 0 { + matrix = matrix.dot(&pauli_matrix); + } else { + matrix = kron(&matrix, &pauli_matrix); + } + } + // now check that this matrix squares to the identity as otherwise the exact value will be + // incorrect. + let matrix_squared = matrix.dot(&matrix); + let diff = &matrix_squared - Array2::::eye(dim); + assert!(diff.opnorm_one().unwrap() < 10. * (dim as f64) * f64::EPSILON); + + let theta = 1. * std::f64::consts::PI * rng.gen::(); + let scaled_matrix = matrix.clone() * c64::new(0., theta); + let expm_computed = expm(&scaled_matrix).unwrap(); + let expm_expected = Array2::::eye(dim) * theta.cos() + c64::new(0., theta.sin()) * matrix; + let comp_diff = &expm_expected - &expm_computed; + + let error = comp_diff.opnorm_one().unwrap() / f64::EPSILON; + assert!(error <= 10.); +} + +/// Test dense matrix exponentiation from random matrices. This works by constructing a random +/// unitary matrix as the eigenvectors and then random eigenvalues. We can then use the f64 +/// exponentiation routine to compute the exponential of the diagonal eigenvalue matrix and then +/// multiply by the eigenvalues to compute the exponentiated matrix. This exact value is compared +/// to the expm calculation by looking at the average error per matrix entry. This test reveals an +/// error that scales with the dimension worse than competitors on the author's laptop (Mac M1 with +/// the Accelerate BLAS/LAPACK backend as of December 10th, 2022). +#[test] +fn test_low_dimension_random_dense_matrix() { + let mut rng = rand::thread_rng(); + let dimensions = 100; + let samps = 500; + let scale = 1.; + let mut results = Vec::new(); + let mut avg_entry_error = Vec::new(); + // Used to control what pade approximation is most likely to be used. + // the smaller the norm the lower the degree used. + for _ in 0..samps { + // Sample a completely random matrix. + let mut matrix: Array2 = Array2::::ones((dimensions, dimensions).f()); + matrix.mapv_inplace(|_| c64::new(rng.gen::() * 1., rng.gen::() * 1.)); + + // Make m positive semidefinite so it has orthonormal eigenvecs. + matrix = matrix.dot(&matrix.t().map(|x| x.conj())); + let (mut eigs, vecs) = matrix.eig().unwrap(); + let adjoint_vecs = vecs.t().clone().mapv(|x| x.conj()); + + // Generate new random eigenvalues (complex, previously m had real eigenvals) + // and a new matrix m + eigs.mapv_inplace(|_| scale * c64::new(rng.gen::(), rng.gen::())); + let new_matrix = vecs.dot(&Array2::from_diag(&eigs)).dot(&adjoint_vecs); + + // compute the exponentiated matrix by exponentiating the eigenvalues + // and doing V e^Lambda V^\dagger + eigs.mapv_inplace(|x| x.exp()); + let eigen_expm = vecs.dot(&Array2::from_diag(&eigs)).dot(&adjoint_vecs); + + // Compute the expm routine, compute error metrics for this sample + let expm_comp = expm(&new_matrix).unwrap(); + let diff = &expm_comp - &eigen_expm; + avg_entry_error.push({ + let tot = diff.map(|x| x.abs()).into_iter().sum::(); + tot / (dimensions * dimensions) as f64 + }); + results.push(diff.opnorm_one().unwrap()); + } + + // compute averages + let avg: f64 = results.iter().sum::() / results.len() as f64; + let avg_entry_diff = avg_entry_error.iter().sum::() / avg_entry_error.len() as f64; + let _std: f64 = f64::powf( + results.iter().map(|x| f64::powi(x - avg, 2)).sum::() / (results.len() - 1) as f64, + 0.5, + ); + + // This may fail at higher dimensions + assert!(avg_entry_diff / f64::EPSILON <= 1000.) +}