# Performance Metrics

The following scripts were used to check errors and runtime.

## Error Calculations

Since parameter sets are identified up to permutation, the following finds the permutation that best reduces the error of the mixing coefficients and applies it to the means and covariances to compute their respective errors.

function computeError(w, true_means, true_covariances, mixing_coefficients, means, covariances, diagonal)
k = size(w)[1]
d = size(true_means)[2]

basis = 1:k

minimum_weight_error = (norm(w - mixing_coefficients), basis)
for i in 2:factorial(k)
permutation = nthperm(basis, i)
mixed_weights = Array{Float64}(undef, size(w))
for j in 1:k
mixed_weights[j, 1:end] = mixing_coefficients[permutation[j], 1:end]
end
weight_error = norm(mixed_weights - w)
(weight_error < minimum_weight_error[1]) && (minimum_weight_error = (weight_error, permutation))
end

permutation = minimum_weight_error[2]

final_mixing_coefficients = Array{Float64}(undef, size(mixing_coefficients))
final_means = Array{Float64}(undef, size(means))

if diagonal
final_covariances = []
for j in 1:k
final_mixing_coefficients[j] = mixing_coefficients[permutation[j]]
final_means[j, 1:end] = means[permutation[j], 1:end]
push!(final_covariances, covariances[permutation[j]][1:end, 1:end])
end
else
final_covariances = Array{Union{Float64}, 3}(undef, (k,d,d))
for j in 1:k
final_mixing_coefficients[j] = mixing_coefficients[permutation[j]]
final_means[j, 1:end] = means[permutation[j], 1:end]
final_covariances[j, 1:end, 1:end] = covariances[permutation[j], 1:end, 1:end]
end
end

mixing_error = norm(final_mixing_coefficients - w)
means_error = norm(final_means - true_means)
covariance_error = norm(final_covariances - true_covariances)

return (mixing_error, means_error, covariance_error)
end

This is similar, but for Algorithm 2 from Estimating Gaussian Mixtures Using Sparse Polynomial Moment Systems all but the means are known so it merely permutes the means to find the minimum error permutation and then returns that error.

function computeErrorAlg2(true_mean, est_mean)
k = size(true_mean)[1]
d = size(true_mean)[2]

basis = 1:k

minimum_error = norm(true_mean - est_mean)
for i in 2:factorial(k)
permutation = nthperm(basis, i)
mixed_means = Array{Float64}(undef, size(true_mean))
for j in 1:k
mixed_means[j, 1:end] = est_mean[permutation[j], 1:end]
end
mean_error = norm(mixed_means - true_mean)
(mean_error < minimum_error) && (minimum_error = mean_error)
end

return minimum_error
end

## Tests

The following tests performance while adding specified levels of noise to the denoised moments. This relies on the above listed computeError function.

function testNoise(d, k, diagonal, noise_levels, reps, known_mixing_coeffs)
passes = []
mix_errs = []
mean_errs = []
covar_errs = []
timings = []

for noise in noise_levels
passing = 0
mixing_error = 0
means_error = 0
covariance_error = 0
times = 0

for i in 1:reps
w, true_means, true_covariances = generateGaussians(d, k, diagonal)
if diagonal
true_first, true_diag = diagonalPerfectMoments(d, k, w, true_means, true_covariances)
num_moments = 3*k + (d-1)*(2*k+1)
else
true_first, true_diag, true_others = densePerfectMoments(d, k, w, true_means, true_covariances)
num_moments = 3*k + (d-1)*(2*k+1) + k*(d^2-d)/2
end

if noise>0
randomness = randn(Int64(num_moments))
randomness = noise/norm(randomness) * randomness

noisy_first = true_first + [0, randomness[1:3*k]...]
noisy_diag = true_diag + reshape(randomness[3*k+1:3*k + (d-1)*(2*k+1)], (d-1,2*k+1))

if !diagonal
noisy_others = Dict{Vector{Int64}, Float64}()
orig = []
new = []
counter = 3*k + (d-1)*(2*k+1) + 1
for (key, moment) in true_others
push!(orig, moment)
push!(new, moment + randomness[counter])
noisy_others[key] = moment + randomness[counter]
counter += 1
end
end
else
noisy_first = true_first
noisy_diag = true_diag
if !diagonal
noisy_others = true_others
end
end

pass = false
if diagonal
if known_mixing_coeffs
timing = @elapsed begin
pass, (mixing_coefficients, means, covariances) = estimate_parameters(d, k, w, noisy_first, noisy_diag)
end
else
timing = @elapsed begin
pass, (mixing_coefficients, means, covariances) = estimate_parameters(d, k, noisy_first, noisy_diag)
end
end
else
if known_mixing_coeffs
timing = @elapsed begin
pass, (mixing_coefficients, means, covariances) = estimate_parameters(d, k, w, noisy_first, noisy_diag, noisy_others)
end
else
timing = @elapsed begin
pass, (mixing_coefficients, means, covariances) = estimate_parameters(d, k, noisy_first, noisy_diag, noisy_others)
end
end
end

if pass == true
passing += 1
(mix, mean, cov) = computeError(w, true_means, true_covariances, mixing_coefficients, means, covariances, diagonal)
mixing_error += mix
means_error += mean
covariance_error += cov
times += timing
end
end
push!(passes, passing)
if passing > 0
push!(mix_errs, mixing_error/passing)
push!(mean_errs, means_error/passing)
push!(covar_errs, covariance_error/passing)
push!(timings, times/passing)
else
push!(mix_errs, nothing)
push!(mean_errs, nothing)
push!(covar_errs, nothing)
push!(timings, nothing)
end
open("noise_test_d" * string(d) * "_k" * string(k) * "_diag" * string(diagonal) * "_knoww" * string(known_mixing_coeffs) * ".txt", "w") do file
write(file, "reps: " * string(reps) * "\nnoise levels: " * string(noise_levels) * "\npasses: " * string(passes) * "\naverage mixing coefficient error: " * string(mix_errs) * "\naverage mean error: " * string(mean_errs) * "\naverage covariance error: " * string(covar_errs) * "\naverage time: " * string(timings) * "\n")
end
println("check")
end
end

Example of use:

k = 3
d = 5
diag = true
known = false
testNoise(d, k, diag, [.1, .01, .001, .0001, 0], 10, known)

The following tests performance with 100 samples per parameter, cycling over dimensions to find candidate mixing coefficients if the solver is unsuccessful with prior mixing coefficient candidates. This relies on the above listed computeError function.

function testSampleAlg2(d, k, reps)
num_params = k*d
num_samples = 200*num_params

passing = 0
means_error = 0
times = 0

w = (1/k)*ones(k)

for i in 1:reps
shared_cov = makeCovarianceMatrix(d, false)
true_mean = randn(k,d)

variances = Array{Float64, 3}(undef, (k,d,d))
for i in 1:k
variances[i, 1:end, 1:end] = copy(shared_cov)
end

sample = getSample(num_samples, w, true_mean, variances)
variances = []

m1, m2 = equalMixCovarianceKnown_moments(k, sample)

timing = @elapsed begin
pass, est_mean = estimate_parameters(k, shared_cov, m1, m2)
end
if pass == true
passing += 1
means_error += computeErrorAlg2(true_mean, est_mean)
times += timing
end
end

if passing > 0
means_error = means_error/passing
times = times/passing
else
means_error = nothing
times = nothing
end

open("alg2_sample_test_d" * string(d) * "_k" * string(k) * ".txt", "w") do file
write(file, "reps: " * string(reps) * "\npasses: " * string(passing) * "\naverage mean error: " * string(means_error) * "\naverage time: " * string(times) * "\n")
end
end

Example of use:

k = 2
diagonal = true
d = 10
reps = 100
testSampleCycle(d, k, diagonal, reps)

The following tests performance of Algorithm 2 from Estimating Gaussian Mixtures Using Sparse Polynomial Moment Systems for denoised moments. This relies on the above listed computeErrorAlg2 function.

function testPerfectAlg2(d, k, reps)
passing = 0
means_error = 0
times = 0

for i in 1:reps
shared_cov = makeCovarianceMatrix(d, false)
true_mean = randn(k,d)
m1, m2 = equalMixCovarianceKnown_moments(k, true_mean, shared_cov)

timing = @elapsed begin
pass, est_mean = estimate_parameters(k, shared_cov, m1, m2)
end
if pass == true
passing += 1
means_error += computeErrorAlg2(true_mean, est_mean)
times += timing
end
end

if passing > 0
means_error = means_error/passing
times = times/passing
else
means_error = nothing
times = nothing
end

open("alg2_perfect_test_d" * string(d) * "_k" * string(k) * ".txt", "w") do file
write(file, "reps: " * string(reps) * "\npasses: " * string(passing) * "\naverage mean error: " * string(means_error) * "\naverage time: " * string(times) * "\n")
end
end

The following tests performance of Algorithm 2 from Estimating Gaussian Mixtures Using Sparse Polynomial Moment Systems for sample moments with 200 samples per parameter. This relies on the above listed computeErrorAlg2 function.

function testSampleAlg2(d, k, reps)
num_params = k*d
num_samples = 200*num_params

passing = 0
means_error = 0
times = 0

w = (1/k)*ones(k)

for i in 1:reps
shared_cov = makeCovarianceMatrix(d, false)
true_mean = randn(k,d)

variances = Array{Float64, 3}(undef, (k,d,d))
for i in 1:k
variances[i, 1:end, 1:end] = copy(shared_cov)
end

sample = getSample(num_samples, w, true_mean, variances)
variances = []

m1, m2 = equalMixCovarianceKnown_moments(k, sample)

timing = @elapsed begin
pass, est_mean = estimate_parameters(k, shared_cov, m1, m2)
end
if pass == true
passing += 1
means_error += computeErrorAlg2(true_mean, est_mean)
times += timing
end
end

if passing > 0
means_error = means_error/passing
times = times/passing
else
means_error = nothing
times = nothing
end

open("alg2_sample_test_d" * string(d) * "_k" * string(k) * ".txt", "w") do file
write(file, "reps: " * string(reps) * "\npasses: " * string(passing) * "\naverage mean error: " * string(means_error) * "\naverage time: " * string(times) * "\n")
end
end