A static class to compute the trace of implicit matrix functions using stochastic Lanczos quadrature method. This class acts as a templated namespace, where the member methods is public and static. The internal private member functions are also static. More...

#include <convergence_tools.h>

Static Public Member Functions
static FlagType	check_convergence (DataType *samples, const IndexType min_num_samples, const IndexType num_inquiries, const IndexType processed_samples_indices, const IndexType num_processed_samples, const DataType confidence_level, const DataType error_atol, const DataType error_rtol, DataType error, IndexType num_samples_used, FlagType *converged)
	Checks if the standard deviation of the set of the cumulative averages of trace estimators converged below the given tolerance.

static void	average_estimates (const DataType confidence_level, const DataType outlier_significance_level, const IndexType num_inquiries, const IndexType max_num_samples, const IndexType num_samples_used, const IndexType processed_samples_indices, DataType *samples, IndexType num_outliers, DataType trace, DataType error)
	Averages the estimates of trace. Removes outliers and reevaluates the error to take into account for the removal of the outliers.

Detailed Description

template<typename DataType>
class ConvergenceTools< DataType >

A static class to compute the trace of implicit matrix functions using stochastic Lanczos quadrature method. This class acts as a templated namespace, where the member methods is public and static. The internal private member functions are also static.

See also: Diagonalization

Definition at line 36 of file convergence_tools.h.

Member Function Documentation

◆ average_estimates()

template<typename DataType >

void ConvergenceTools< DataType >::average_estimates	(	const DataType	confidence_level,
		const DataType	outlier_significance_level,
		const IndexType	num_inquiries,
		const IndexType	max_num_samples,
		const IndexType *	num_samples_used,
		const IndexType *	processed_samples_indices,
		DataType **	samples,
		IndexType *	num_outliers,
		DataType *	trace,
		DataType *	error
	)

static

Averages the estimates of trace. Removes outliers and reevaluates the error to take into account for the removal of the outliers.

Note: The elimination of outliers does not affect the elements of samples array, rather it only affects the re-evaluation of trace and error arrays.

Parameters

[in]	confidence_level	The confidence level of the error, which is a number between `0` and `1`. This affects the scale of `error`.
[in]	outlier_significance_level	One minus the confidence level of the uncertainty band of the outlier. This is a number between `0` and `1`. Confidence level of outleir and significance level of outlier are commlement of each other.
[in]	num_inquiries	The number of batches of parameters. This function computes `num_inquiries` values of trace corresponding to different batch of parameters of the linear operator `A`. Hence, the number of output trace is `num_inquiries`. Hence, it is the number of columns of the output array `samples`.
[in]	max_num_samples	The number of times that the trace estimation is repeated. The output trace value is the average of the samples. Hence, this is the number of rows of the output array `samples`. Larger number of samples leads to a better trace estimation. The computational const linearly increases with number of samples.
[in]	num_samples_used	1D array of the size of the number of columns of `samples`. Each element indicates how many iterations were used till convergence is reached for each column of the `samples`. The number of iterations should be a number between `min_num_samples` and `max_num_samples`.
[in]	processed_samples_indices	A 1D array indicating the processing order of rows of the `samples`. In paralleli processing, this order of processing the rows of `samples` is not necessarly sequential.
[out]	samples	2D array of all estimated trace samples. The shape of this array is (max_num_samples*num_inquiries). The average of the rows is also given in `trace` array.
[out]	num_outliers	1D array with the size of number of columns of `samples`. Each element indicates how many rows of the `samples` array were outliers and were removed during averaging rows of `samples`.
[out]	trace	The output trace of size `num_inquiries`. These values are the average of the rows of `samples` array.
[out]	error	The error of estimation of trace, which is the standard deviation of the rows of `samples` array. The size of this array is `num_inquiries`.

Definition at line 289 of file convergence_tools.cpp.

{
    IndexType i;
    IndexType j;
    DataType summand;
    DataType mean;
    DataType mean_abs;
    DataType std;
    DataType mean_discrepancy;
    DataType outlier_half_interval;
 
    // Flag which samples are outliers
    FlagType* outlier_indices = new FlagType[max_num_samples];
 
    // Quantile of normal distribution (usually known as the "z" coefficient)
    DataType error_z_score = std::sqrt(2) * erf_inv(confidence_level);
 
    // Confidence level of outlier is the complement of significance level
    DataType outlier_confidence_level = 1.0 - outlier_significance_level;
 
    // Quantile of normal distribution area where is not considered as outlier
    DataType outlier_z_score = std::sqrt(2.0) * \
        erf_inv(outlier_confidence_level);
 
    for (j=0; j < num_inquiries; ++j)
    {
        // Initialize outlier indices for each column of samples
        for (i=0; i < max_num_samples; ++i)
        {
            outlier_indices[i] = 0;
        }
        num_outliers[j] = 0;
 
        // Compute mean of the j-th column
        summand = 0.0;
        for (i=0; i < num_samples_used[j]; ++i)
        {
            summand += samples[processed_samples_indices[i]][j];
        }
        mean = summand / num_samples_used[j];
 
        // Compute mean of the absolute values of j-th column
        summand = 0.0;
        for (i=0; i < num_samples_used[j]; ++i)
        {
            summand += std::abs(samples[processed_samples_indices[i]][j]);
        }
        mean_abs = summand / num_samples_used[j];
 
        // If mean of absolute values is zero, do not scale data with it as it
        // causes divide by zero. In this case, all sample data are zero, and
        // there is no need to scale them.
        DataType zero = 0.0;
        if (c_arithmetics::is_equal(mean_abs, zero))
        {
            mean_abs = 1.0;
        }
 
        // Compute std of the j-th column
        if (num_samples_used[j] > 1)
        {
            summand = 0.0;
            for (i=0; i < num_samples_used[j]; ++i)
            {
                mean_discrepancy = \
                    samples[processed_samples_indices[i]][j] - mean;
 
                // Normalize to the mean of absolute values to avoid underflow
                // and overflow, but later, re-scale it back. The underflow or
                // overflow is caused by taking the square two lines below.
                mean_discrepancy /= mean_abs;
 
                summand += mean_discrepancy * mean_discrepancy;
            }
            std = std::sqrt(summand / (num_samples_used[j] - 1.0));
 
            // Re-scale back by mean of absolute values
            std *= mean_abs;
        }
        else
        {
            std = INFINITY;
        }
 
        // Outlier half interval
        outlier_half_interval = outlier_z_score * std;
 
        // Find outliers by the difference of each element from mean
        for (i=0; i < num_samples_used[j]; ++i)
        {
            mean_discrepancy = samples[processed_samples_indices[i]][j] - mean;
            if (std::abs(mean_discrepancy) > outlier_half_interval)
            {
                // Outlier detected
                outlier_indices[i] = 1;
                num_outliers[j] += 1;
            }
        }
 
        // Re-evaluate mean but leave out outliers
        summand = 0.0;
        for (i=0; i < num_samples_used[j]; ++i)
        {
            if (outlier_indices[i] == 0)
            {
                summand += samples[processed_samples_indices[i]][j];
            }
        }
        mean = summand / (num_samples_used[j] - num_outliers[j]);
 
        // Re-evaluate std but leave out outliers
        if (num_samples_used[j] > 1 + num_outliers[j])
        {
            summand = 0.0;
            for (i=0; i < num_samples_used[j]; ++i)
            {
                if (outlier_indices[i] == 0)
                {
                    mean_discrepancy = \
                        samples[processed_samples_indices[i]][j] - mean;
 
                // Normalize to the mean of absolute values to avoid underflow
                // and overflow, but later, re-scale it back. The underflow or
                // overflow is caused by taking the square two lines below.
                mean_discrepancy /= mean_abs;
 
                    summand += mean_discrepancy * mean_discrepancy;
                }
            }
            std = std::sqrt(
                summand/(num_samples_used[j] - num_outliers[j] - 1.0));
 
            // Re-scale back by mean of absolute values
            std *= mean_abs;
        }
        else
        {
            std = INFINITY;
        }
 
        // trace and its error
        trace[j] = mean;
        error[j] = error_z_score * std / \
            std::sqrt(num_samples_used[j] - num_outliers[j]);
    }
 
    delete[] outlier_indices;
}

References erf_inv(), and c_arithmetics::is_equal().

Referenced by cTraceEstimator< DataType >::c_trace_estimator(), and cuTraceEstimator< DataType >::cu_trace_estimator().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ check_convergence()

template<typename DataType >

FlagType ConvergenceTools< DataType >::check_convergence	(	DataType **	samples,
		const IndexType	min_num_samples,
		const IndexType	num_inquiries,
		const IndexType *	processed_samples_indices,
		const IndexType	num_processed_samples,
		const DataType	confidence_level,
		const DataType	error_atol,
		const DataType	error_rtol,
		DataType *	error,
		IndexType *	num_samples_used,
		FlagType *	converged
	)

static

Checks if the standard deviation of the set of the cumulative averages of trace estimators converged below the given tolerance.

The convergence criterion for each trace inquiry is if:

            standard_deviation < max(rtol * average[i], atol)

        where \c rtol and \c atol are relative and absolute tolerances,
        respectively. If this criterion is satisfied for *all* trace
        inquiries, this function returns \c 1, otherwise \c 0.

Parameters

[in]	samples	A 2D array of the estimated trace. This array has the shape (max_num_samples,num_inquiries). Each column of this array is the estimated trace values for a different trace inquiry. The rows are Monte-Carlo samples of the estimation of trace. The cumulative average of the rows are expected to converge. Note that because of parallel processing, the rows of this array are not filled sequentially. The list of filled rows are stored in `processed_samples_indices`.
[in]	min_num_samples	Minimum number of sample iterations.
[in]	num_inquiries	Number of columns of `samples` array.
[in]	processed_samples_indices	A 1D array indicating the processing order of rows of the `samples`. In parallel processing, this order of processing the rows of `samples` is not necessarly sequential.
[in]	num_processed_samples	A counter that keeps track of how many samples were processed so far in the iterations.
[in]	confidence_level	The confidence level of the error, which is a number between `0` and `1`. This affects the scale of `error`.
[in]	error_atol	Absolute tolerance criterion to be compared with the standard deviation of the cumulative averages. If `error_atol` is zero, it is ignored, and only `rtol` is used.
[in]	error_rtol	Relative tolerance criterion to be compared with the standard deviation of the cumulative averages. If `error_rtol` is zero, it is ignored, and only `error_atol` is used.
[out]	error	The error of estimation of trace, which is the standard deviation of the rows of `samples` array. The size of this array is `num_inquiries`.
[out]	num_samples_used	1D array of the size of the number of columns of `samples`. Each element indicates how many iterations were used till convergence is reached for each column of the `samples`. The number of iterations should be a number between `min_num_samples` and `max_num_samples`.
[out]	converged	1D array of the size of the number of columns of `samples`. Each element indicates which column of `samples` has converged to the tolerance criteria. Normally, if the `num_samples_used` is less than `max_num_samples`, it indicates that the convergence has reached. the rows of `samples` array. The size of this array is `num_inquiries`.

Returns

A signal to indicate the status of computation:

1 indicates successful convergence within the given tolerances was met. Convergence is achieved when all elements of convergence array are below convergence_atol or convergence_rtol times trace.
0 indicates the convergence criterion was not met for at least one of the trace inquiries.

Definition at line 96 of file convergence_tools.cpp.

{
    FlagType all_converged;
    IndexType j;
 
    // If number of processed samples are not enough, set to not converged yet.
    // This is essential since in the first few iterations, the standard
    // deviation of the cumulative averages are still too small.
    if (num_processed_samples < min_num_samples)
    {
        // Skip computing error. Fill outputs with trivial initial values
        for (j=0; j < num_inquiries; j++)
        {
            error[j] = INFINITY;
            converged[j] = 0;
            num_samples_used[j] = num_processed_samples;
        }
        all_converged = 0;
        return all_converged;
    }
 
    IndexType i;
    DataType summand;
    DataType mean;
    DataType mean_abs;
    DataType std;
    DataType mean_discrepancy;
    DataType data;
 
    // Quantile of normal distribution (usually known as the "z" coefficient)
    DataType standard_z_score = std::sqrt(2) * \
        static_cast<DataType>(erf_inv(static_cast<double>(confidence_level)));
 
    // For each column of samples, compute error of all processed rows
    for (j=0; j < num_inquiries; ++j)
    {
        // Do not check convergence if j-th column already converged
        if (converged[j] == 0)
        {
            // mean of j-th column using all processed rows of j-th column
            summand = 0.0;
            for (i=0; i < num_processed_samples; ++i)
            {
                summand += samples[processed_samples_indices[i]][j];
            }
            mean = summand / num_processed_samples;
 
            // mean of absolute values of j-th column using all processed rows
            // of j-th column
            summand = 0.0;
            for (i=0; i < num_processed_samples; ++i)
            {
                summand += std::abs(samples[processed_samples_indices[i]][j]);
            }
            mean_abs = summand / num_processed_samples;
 
            // If mean of absolute values is zero, do not scale data with it as
            // it causes divide by zero. In this case, all sample data are
            // zero, and there is no need to scale them.
            DataType zero = 0.0;
            if (c_arithmetics::is_equal(mean_abs, zero))
            {
                mean_abs = 1.0;
            }
 
            // std of j-th column using all processed rows of j-th column
            if (num_processed_samples > 1)
            {
                summand = 0.0;
                for (i=0; i < num_processed_samples; ++i)
                {
                    data = samples[processed_samples_indices[i]][j];
                    mean_discrepancy = data - mean;
 
                    // Normalize to the mean of absolute values to avoid
                    // underflow and overflow, but later, re-scale it back. The
                    // underflow or overflow is caused by taking the square two
                    // lines below.
                    mean_discrepancy /= mean_abs;
                    
                    summand += mean_discrepancy * mean_discrepancy;
                }
                std = std::sqrt(summand / (num_processed_samples - 1.0));
 
                // Re-scale back by mean of absolute values
                std *= mean_abs;
            }
            else
            {
                std = INFINITY;
            }
 
            // Compute error based of std and confidence level
            error[j] = standard_z_score * std / \
                std::sqrt(num_processed_samples);
 
            // Check error with atol and rtol to find if j-th column converged
            if (error[j] < std::max(error_atol, error_rtol*mean))
            {
                converged[j] = 1;
            }
 
            // Update how many samples used so far to average j-th column
            num_samples_used[j] = num_processed_samples;
        }
    }
 
    // Check convergence is reached for all columns (all inquiries)
    all_converged = 1;
    for (j=0; j < num_inquiries; ++j)
    {
        if (converged[j] == 0)
        {
            // The j-th column not converged.
            all_converged = 0;
            break;
        }
    }
 
    return all_converged;
}