Utils & Metrics

betainc

pyldl.algorithms.utils.betainc(alpha, beta, x)

Compute the regularized incomplete beta function.

Parameters:

• alpha (np.ndarray or tensor) – First beta-shape parameter.

• beta (np.ndarray or tensor) – Second beta-shape parameter.

• x (np.ndarray or tensor) – Evaluation points.

Returns:

Beta cumulative distribution values.

Return type:

np.ndarray or tensor

binaryzation

pyldl.algorithms.utils.binaryzation(D: ndarray, method='threshold', param: any | None = None) → ndarray

Transform label distribution matrix to logical label matrix.

Parameters:

• D (np.ndarray) – Label distribution matrix (shape: \([n, c]\)).

• method ({'threshold', 'topk'}, optional) –

  Type of binaryzation method, defaults to ‘threshold’. The options are ‘threshold’ and ‘topk’, which can refer to:

  [BIN-KWT+24]

  Zhiqiang Kou, Jing Wang, Jiawei Tang, Yuheng Jia, Boyu Shi, and Xin Geng. Exploiting multi-label correlation in label distribution learning. In Proceedings of the International Joint Conference on Artificial Intelligence, 4326–4334. 2024. URL: https://doi.org/10.24963/ijcai.2024/478.

• param (any, optional) – Parameter of binaryzation method, defaults to None. If None, the default value is .5 for ‘threshold’ and \(\lfloor c / 2 \rfloor\) for ‘topk’.

Returns:

Logical label matrix (shape: \([n, c]\)).

Return type:

np.ndarray

digamma

pyldl.algorithms.utils.digamma(x)

Approximate the digamma function of an input array or tensor.

Parameters:

x (np.ndarray or tensor) – Input array or tensor.

Returns:

Element-wise digamma values.

Return type:

np.ndarray or tensor

estimate_alpha

pyldl.algorithms.utils.estimate_alpha(D, max_iterations: int = 100, convergence_criterion=1e-07)

Estimate Dirichlet concentration parameters from label distributions.

Parameters:

• D (np.ndarray or tensor) – Label distribution matrix.

• max_iterations (int) – Maximum number of fixed-point iterations, defaults to 100.

• convergence_criterion (float) – Convergence threshold, defaults to 1e-7.

Returns:

Estimated Dirichlet concentration vector.

Return type:

np.ndarray or tensor

gammaln

pyldl.algorithms.utils.gammaln(x)

Compute the natural logarithm of the gamma function using a Lanczos approximation.

Parameters:

x (np.ndarray or tensor) – Input array or tensor.

Returns:

Element-wise log-gamma values.

Return type:

np.ndarray or tensor

inv_digamma

pyldl.algorithms.utils.inv_digamma(y, max_iterations: int = 2)

Compute the inverse digamma function using Newton iterations.

Parameters:

• y (np.ndarray or tensor) – Digamma values.

• max_iterations (int) – Maximum number of Newton iterations, defaults to 2.

Returns:

Inverse-digamma values.

Return type:

np.ndarray or tensor

kernel

pyldl.algorithms.utils.kernel(X, Y=None, gamma: float | None = None)

Compute an RBF kernel matrix between rows of two matrices.

Parameters:

• X (np.ndarray) – First input matrix.

• Y (np.ndarray or None) – Second input matrix; if None, use X, defaults to None.

• gamma (float or None) – RBF kernel coefficient; if None, estimate it from pairwise distances, defaults to None.

Returns:

RBF kernel matrix.

Return type:

np.ndarray

kl_divergence\(\downarrow\)

pyldl.algorithms.utils.kl_divergence(D, D_pred)

Kullback-Leibler divergence. It is defined as:

\[\text{KLD}(\boldsymbol{u}, \boldsymbol{v}) = \sum^c_{j=1}u_j \ln \frac{u_j}{v_j}\text{.}\]

locass_proj

pyldl.algorithms.utils.locass_proj(D, k)

non_diagonal

pyldl.algorithms.utils.non_diagonal(X)

Set the diagonal entries of a square matrix to zero.

Parameters:

X (np.ndarray) – Input square matrix.

Returns:

Matrix with zero diagonal.

Return type:

np.ndarray

normalize

pyldl.algorithms.utils.normalize(D: ndarray) → ndarray

Normalize each row of a matrix to sum to one.

Parameters:

D (np.ndarray) – Input matrix.

Returns:

Row-normalized matrix.

Return type:

np.ndarray

pairwise_cosine

pyldl.algorithms.utils.pairwise_cosine(X, Y=None, mode: str = 'similarity')

Compute pairwise cosine similarities or distances between rows of two matrices.

Parameters:

• X (np.ndarray) – First input matrix.

• Y (np.ndarray or None) – Second input matrix; if None, use X, defaults to None.

• mode ({'similarity', 'distance'}) – Return mode, either 'similarity' or 'distance', defaults to 'similarity'.

Returns:

Pairwise cosine similarity or distance matrix.

Return type:

np.ndarray

pairwise_euclidean

pyldl.algorithms.utils.pairwise_euclidean(X, Y=None)

Compute pairwise Euclidean distances between rows of two matrices.

Parameters:

• X (np.ndarray) – First input matrix.

• Y (np.ndarray or None) – Second input matrix; if None, use X, defaults to None.

Returns:

Pairwise Euclidean distance matrix.

Return type:

np.ndarray

pairwise_pearsonr

pyldl.algorithms.utils.pairwise_pearsonr(X, Y=None)

Compute pairwise Pearson correlation coefficients between rows of two matrices.

Parameters:

• X (np.ndarray) – First input matrix.

• Y (np.ndarray or None) – Second input matrix; if None, use X, defaults to None.

Returns:

Pairwise Pearson correlation matrix.

Return type:

np.ndarray

proj

pyldl.algorithms.utils.proj(D: ndarray) → ndarray

This approach is proposed in paper [Con16].

Parameters:

D (np.ndarray) – Input matrix.

Returns:

The projection onto the probability simplex.

Return type:

np.ndarray

shannon_entropy

pyldl.algorithms.utils.shannon_entropy(D: ndarray)

Compute the Shannon entropy of a label distribution matrix.

Parameters:

D (np.ndarray) – Label distribution matrix.

Returns:

Shannon entropy.

Return type:

float

soft_thresholding

pyldl.algorithms.utils.soft_thresholding(A: ndarray, tau: float) → ndarray

Soft thresholding operation. It is defined as \(\text{soft}(\boldsymbol{A}, \tau) = \text{sgn}(\boldsymbol{A}) \odot \max\lbrace \lvert \boldsymbol{A} \rvert - \tau, 0 \rbrace\), where \(\odot\) denotes element-wise multiplication.

Parameters:

• A (np.ndarray) – Matrix \(\boldsymbol{A}\).

• tau (float) – \(\tau\).

Returns:

The result of soft thresholding operation.

Return type:

np.ndarray

softmax

pyldl.algorithms.utils.softmax(D: ndarray) → ndarray

Apply the softmax function to each row of a matrix.

Parameters:

D (np.ndarray) – Input matrix.

Returns:

Row-wise softmax values.

Return type:

np.ndarray

solvel21

pyldl.algorithms.utils.solvel21(A: ndarray, tau: float) → ndarray

This approach is proposed in paper [CY14].

The solution to the optimization problem \(\mathop{\arg\min}_{\boldsymbol{X}} \Vert \boldsymbol{X} - \boldsymbol{A} \Vert_\text{F}^2 + \tau \Vert \boldsymbol{X} \Vert_{2,1}\) is given by the following formula:

\[\begin{split}\vec{x}_{\bullet j}^{\ast} = \left\{ \begin{aligned} & \frac{\Vert \vec{a}_{\bullet j} \Vert - \tau}{\Vert \vec{a}_{\bullet j} \Vert} \vec{a}_{\bullet j}, & \tau \le \Vert \vec{a}_{\bullet j} \Vert \\ & 0, & \text{otherwise} \end{aligned} \right.\text{.}\end{split}\]

where \(\vec{x}_{\bullet j}\) is the \(j\)-th column of matrix \(\boldsymbol{X}\), and \(\vec{a}_{\bullet j}\) is the \(j\)-th column of matrix \(\boldsymbol{A}\).

Parameters:

• A (np.ndarray) – Matrix \(\boldsymbol{A}\).

• tau (float) – \(\tau\).

Returns:

The solution to the optimization problem.

Return type:

np.ndarray

svt

pyldl.algorithms.utils.svt(A: ndarray, tau: float) → ndarray

Singular value thresholding (SVT) is proposed in paper [CCS10].

The solution to the optimization problem \(\mathop{\arg\min}_{\boldsymbol{X}} \Vert \boldsymbol{X} - \boldsymbol{A} \Vert_\text{F}^2 + \tau \Vert \boldsymbol{X} \Vert_{\ast}\) is given by \(\boldsymbol{U} \max \lbrace \boldsymbol{\Sigma} - \tau, 0 \rbrace \boldsymbol{V}^\top\), where \(\boldsymbol{A} = \boldsymbol{U} \boldsymbol{\Sigma} \boldsymbol{V}^\top\) is the singular value decomposition of matrix \(\boldsymbol{A}\).

Parameters:

• A (np.ndarray) – Matrix \(\boldsymbol{A}\).

• tau (float) – \(\tau\).

Returns:

The solution to the optimization problem.

Return type:

np.ndarray

trigamma

pyldl.algorithms.utils.trigamma(x)

Compute the trigamma function.

Parameters:

x (tensor) – Input tensor.

Returns:

Element-wise trigamma values.

Return type:

tensor

artificial

pyldl.utils.artificial(X, a=1.0, b=0.5, c=0.2, d=1.0, w1=array([[4., 2., 1.]]), w2=array([[1., 2., 4.]]), w3=array([[1., 4., 2.]]), lambda1=0.01, lambda2=0.01)

Generate artificial label distributions from polynomial feature interactions. The generation process is provided in [LDL-Gen16].

Parameters:

• X (np.ndarray) – Feature matrix.

• a – Coefficient of the linear term, defaults to 1.

• b – Coefficient of the quadratic term, defaults to 0.5.

• c – Coefficient of the cubic term, defaults to 0.2.

• d – Constant term, defaults to 1.

• w1 – Weights for the first label, defaults to [[4., 2., 1.]].

• w2 – Weights for the second label, defaults to [[1., 2., 4.]].

• w3 – Weights for the third label, defaults to [[1., 4., 2.]].

• lambda1 – Coupling coefficient for the second label, defaults to 0.01.

• lambda2 – Coupling coefficient for the third label, defaults to 0.01.

Returns:

Artificial label distribution matrix.

Return type:

np.ndarray

download_dataset

pyldl.utils.download_dataset(name, dataset_path)

Download a named .mat dataset from the PyLDL GitHub repository.

Parameters:

• name (str) – Name of the dataset (without the .mat extension).

• dataset_path (str) – Local path where the downloaded dataset is saved.

emphasize

pyldl.utils.emphasize(D, rate=0.5, **kwargs)

gaussian_noise

pyldl.utils.gaussian_noise(D: ndarray, mode=None, mean: float = 0.0, std: float = 0.1, p: float = 5.0)

Add Gaussian noise to label distributions and renormalize them.

Parameters:

• D (np.ndarray) – Label distribution matrix.

• mode ({None, 'gcia'}) –

  Noise mode, defaults to None. Mode ‘gcia’ refers to Generative Calibration of Inaccurate Annotators. See the following reference for details:

  [NOISE-HLLJ24]

  Liang He, Yunan Lu, Weiwei Li, and Xiuyi Jia. Generative calibration of inaccurate annotation for label distribution learning. In Proceedings of the AAAI Conference on Artificial Intelligence, 12394–12401. 2024. URL: https://doi.org/10.1609/aaai.v38i11.29131.

• mean (float) – Mean of the noise in the default mode, defaults to 0.

• std (float) – Standard deviation of the noise in the default mode, defaults to 0.1.

• p (float) – Annotator professionalism level in GCIA mode, defaults to 5.; larger values produce more noise.

Returns:

Noisy label distribution matrix.

Return type:

np.ndarray

load_dataset

pyldl.utils.load_dataset(name, dir='datasets', mode='D')

Load a dataset from a local .mat file, downloading it when necessary. If the dataset is not found in the directory, and it is available in the PyLDL GitHub repository, it will be downloaded automatically.

Parameters:

• name (str) – Name of the dataset (without the .mat extension).

• dir (str) – Directory where the dataset is stored (or will be downloaded to), defaults to ‘datasets’.

• mode ({'D', 'G'}) – Mode of the dataset to load, defaults to ‘D’. ‘D’ is for label distribution, and ‘G’ is for generalized label distribution.

make_ldl

pyldl.utils.make_ldl(n_samples=200, random_state=None, **kwargs)

Generate random features and artificial label distributions using the artificial() function.

Parameters:

• n_samples (int) – Number of samples to generate, defaults to 200.

• random_state (int or None) – Seed for NumPy’s random number generator, defaults to None.

• kwargs – Additional keyword arguments passed to artificial().

Returns:

Generated feature matrix and label distribution matrix.

Return type:

tuple[np.ndarray, np.ndarray]

plot_artificial

pyldl.utils.plot_artificial(grid_size=50, model=None, file_name=None, *, noise=False, noise_func_args=None, **kwargs)

Plot predictions or perturbations on the artificial LDL dataset as a 3D color surface.

Parameters:

• grid_size (int) – Number of points along each feature axis, defaults to 50.

• model (BaseLDL or BaseLE or None) – Optional LDL or LE model used to generate predictions, defaults to None.

• file_name (str or None) – Output file name without extension; if None, display the figure, defaults to None.

• noise (bool) – Whether to apply the default noise transformations, defaults to False.

• noise_func_args (list or None) – Noise functions and their keyword arguments, defaults to None.

• kwargs – Additional keyword arguments passed to artificial().

random_exchange

pyldl.utils.random_exchange(D, rate=0.2, weighted=False, return_mask=True)

random_missing

pyldl.utils.random_missing(D, rate=0.8, weighted=False, return_mask=True)

Randomly set entries of label distributions to zero. This setting is useful for simulating incomplete label distributions. See [INCOM-LDL-XZ17] for details.

Parameters:

• D (np.ndarray) – Label distribution matrix.

• rate (float) – Proportion of entries to remove, in the range (0, 1), defaults to 0.8.

• weighted (bool) –

  Whether to select entries according to value-based weights, defaults to False. This setting can refer to:

  [NOISE-LC24]

  Xiang Li and Songcan Chen. No regularization is needed: efficient and effective incomplete label distribution learning. In Proceedings of the International Joint Conference on Artificial Intelligence, 4470–4478. 2024. URL: https://doi.org/10.24963/ijcai.2024/494.

• return_mask (bool) – Whether to return the missing-entry mask, defaults to True.

Returns:

The incomplete distribution, optionally together with its boolean mask.

Return type:

np.ndarray or tuple[np.ndarray, np.ndarray]

regressor2ldl

pyldl.utils.regressor2ldl(regressor) → BaseLDL

Wrap a scikit-learn regressor as a PyLDL label distribution learner.

Parameters:

regressor (sklearn.base.RegressorMixin) – Scikit-learn regressor used for multi-output prediction.

Returns:

A PyLDL model backed by the supplied regressor.

Return type:

BaseLDL

accuracy\(\uparrow\)

pyldl.metrics.accuracy(y, y_pred)

canberra\(\downarrow\)

pyldl.metrics.canberra(D, D_pred)

Canberra distance. It is defined as:

\[\text{Can.}(\boldsymbol{u}, \boldsymbol{v}) = \sum^c_{j=1}\frac{\left\vert u_j - v_j \right\vert}{u_j + v_j}\text{.}\]

chebyshev\(\downarrow\)

pyldl.metrics.chebyshev(D, D_pred)

Chebyshev distance. It is defined as:

\[\text{Cheby.}(\boldsymbol{u}, \boldsymbol{v}) = \max_j \left\vert u_j - v_j \right\vert\text{.}\]

chi2\(\downarrow\)

pyldl.metrics.chi2(D, D_pred)

Chi-squared distance. It is defined as:

\[\chi^2(\boldsymbol{u}, \boldsymbol{v}) = \sum^c_{j=1}\frac{\left( u_j - v_j \right)^2}{u_j + v_j}\text{.}\]

clark\(\downarrow\)

pyldl.metrics.clark(D, D_pred)

Clark distance. It is defined as:

\[\text{Clark}(\boldsymbol{u}, \boldsymbol{v}) = \sqrt{\sum^c_{j=1}\frac{\left( u_j - v_j \right)^2}{\left( u_j + v_j \right)^2}}\text{.}\]

cosine\(\uparrow\)

pyldl.metrics.cosine(D, D_pred)

Cosine similarity. It is defined as:

\[\text{Cosine}(\boldsymbol{u}, \boldsymbol{v}) = \frac{\sum^c_{j=1}u_j v_j}{\sqrt{\sum^c_{j=1}u_j^2}\sqrt{\sum^c_{j=1}v_j^2}}\text{.}\]

divisiveness_error\(\downarrow\)

pyldl.metrics.divisiveness_error(D, D_pred, pos, neg)

dpa\(\uparrow\)

pyldl.metrics.dpa(D, D_pred)

Degree percentile average (DPA) is proposed in paper [LDL-JQLL24]. It is defined as:

\[\text{DPA}(\boldsymbol{u}, \boldsymbol{v}) = \frac{1}{c} \sum_{j=1}^{c} u_j \rho(v_j)\text{,}\]

where \(\rho(\cdot)\) is the rank of the element in the vector.

error_probability\(\downarrow\)

pyldl.metrics.error_probability(D, D_pred)

Error probability. It is defined as:

\[\text{Err. prob.}(\boldsymbol{u}, \boldsymbol{v}) = 1 - u_{\arg\max(\boldsymbol{v})}\text{.}\]

euclidean\(\downarrow\)

pyldl.metrics.euclidean(D, D_pred)

Euclidean distance. It is defined as:

\[\text{Eucl.}(\boldsymbol{u}, \boldsymbol{v}) = \sqrt{\sum^c_{j=1}\left( u_j - v_j \right)^2}\text{.}\]

fidelity\(\uparrow\)

pyldl.metrics.fidelity(D, D_pred)

Fidelity similarity. It is defined as:

\[\text{Fid.}(\boldsymbol{u}, \boldsymbol{v}) = \sum^c_{j=1} \sqrt{u_j v_j}\text{.}\]

hamming\(\downarrow\)

pyldl.metrics.hamming(L, L_pred)

intersection\(\uparrow\)

pyldl.metrics.intersection(D, D_pred)

Intersection similarity. It is defined as:

\[\text{Int.}(\boldsymbol{u}, \boldsymbol{v}) = \sum^c_{j=1} \min\left(u_j, v_j\right)\text{.}\]

jaccard\(\uparrow\)

pyldl.metrics.jaccard(L, L_pred)

js_divergence\(\downarrow\)

pyldl.metrics.js_divergence(D, D_pred)

Jensen-Shannon divergence. It is defined as:

\[\text{JSD}(\boldsymbol{u}, \boldsymbol{v}) = \frac{1}{2}\text{KLD}\left(\boldsymbol{u} \bigg\Vert \frac{1}{2}(\boldsymbol{u} + \boldsymbol{v}) \right) + \frac{1}{2}\text{KLD}\left(\boldsymbol{v} \bigg\Vert \frac{1}{2}(\boldsymbol{u} + \boldsymbol{v}) \right)\text{.}\]

jsd_dirichlet_expectation

pyldl.metrics.jsd_dirichlet_expectation(alpha, beta)

Estimate the Jensen-Shannon divergence expectation of two Dirichlet distributions.

jsd_dirichlet_variance

pyldl.metrics.jsd_dirichlet_variance(alpha, beta)

Estimate the Jensen-Shannon divergence variance of two Dirichlet distributions.

kendall\(\uparrow\)

pyldl.metrics.kendall(D, D_pred, transpose=False)

Kendall’s rank correlation coefficient. It is defined as:

\[\text{Ken.}(\boldsymbol{u}, \boldsymbol{v}) = \frac{2 \sum_{j < k} \text{sgn}(u_j - u_k) \text{sgn}(v_j - v_k) }{c (c-1)}\text{.}\]

kendallT\(\uparrow\)

pyldl.metrics.kendallT(G, G_pred)

match_m\(\uparrow\)

pyldl.metrics.match_m(D, D_pred, m=None)

max_roc_auc\(\uparrow\)

pyldl.metrics.max_roc_auc(D, D_pred)

mean_absolute_error\(\downarrow\)

pyldl.metrics.mean_absolute_error(D, D_pred, mode='macro')

mean_squared_error\(\downarrow\)

pyldl.metrics.mean_squared_error(D, D_pred, mode='macro')

mu\(\uparrow\)

pyldl.metrics.mu(D, D_pred)

The \(\mu\) metric is proposed in paper [LDL-LWLJ25]. Its KL-divergence-based form is defined as:

\[\mu(\boldsymbol{U}, \boldsymbol{V}) = \frac{1}{\delta_0} \int_0^{\delta_0} \frac{1}{n} \sum_{i=1}^{n} \mathbb{I} (\text{KLD}(\boldsymbol{u}_i, \boldsymbol{v}_i) \le \delta) \mathrm{d}\delta\text{,}\]

where \(\delta_0 = \mathbb{E}_n[\text{KLD}(\boldsymbol{u}_i, \boldsymbol{c})]\) and \(\boldsymbol{c}\) is a uniform vector.

nu1\(\uparrow\)

pyldl.metrics.nu1(D, D_pred, *args, **kwargs)

nu2\(\uparrow\)

pyldl.metrics.nu2(D, D_pred, *args, **kwargs)

nu3\(\uparrow\)

pyldl.metrics.nu3(D, D_pred)

ood_error\(\downarrow\)

pyldl.metrics.ood_error(G, G_pred)

precision\(\uparrow\)

pyldl.metrics.precision(y, y_pred)

score

pyldl.metrics.score(target: ndarray, pred: ndarray, metrics: list | None = None, return_dict: bool = False)

Evaluate one or more metrics for targets and predictions.

sensitivity\(\uparrow\)

pyldl.metrics.sensitivity(y, y_pred)

sorensen\(\downarrow\)

pyldl.metrics.sorensen(D, D_pred)

Sørensen's distance. It is defined as:

\[\text{S}\phi\text{ren.}(\boldsymbol{u}, \boldsymbol{v}) = \frac{\sum^c_{j=1}\left\vert u_j - v_j \right\vert}{\sum^c_{j=1}\left( u_j + v_j \right)}\text{.}\]

spearman\(\uparrow\)

pyldl.metrics.spearman(D, D_pred, transpose=False)

Spearman’s rank correlation coefficient. It is defined as:

\[\text{Spear.}(\boldsymbol{u}, \boldsymbol{v}) = 1 - \frac{6 \sum_{j=1}^{c} (\rho(u_j) - \rho(v_j))^2 }{c(c^2 - 1)}\text{,}\]

where \(\rho(\cdot)\) is the rank of the element in the vector.

spearmanT\(\uparrow\)

pyldl.metrics.spearmanT(G, G_pred)

specificity\(\uparrow\)

pyldl.metrics.specificity(y, y_pred)

subset_accuracy\(\uparrow\)

pyldl.metrics.subset_accuracy(L, L_pred)

top_k\(\uparrow\)

pyldl.metrics.top_k(D, D_pred, k=None, mode='f1_score')

wave_hedges\(\downarrow\)

pyldl.metrics.wave_hedges(D, D_pred)

Wave-Hedges distance. It is defined as:

\[\text{WHD}(\boldsymbol{u}, \boldsymbol{v}) = \sum^c_{j=1}\frac{\left| u_j - v_j \right|}{\max (u_j, v_j)}\text{.}\]

worst_kl_divergence

pyldl.metrics.worst_kl_divergence(D: ndarray)

youden_index\(\uparrow\)

pyldl.metrics.youden_index(y, y_pred)

zero_one_loss\(\downarrow\)

pyldl.metrics.zero_one_loss(D, D_pred)

0/1 loss. It is defined as:

\[\text{0/1 loss}(\boldsymbol{u}, \boldsymbol{v}) = \delta(\arg\max(\boldsymbol{u}), \arg\max(\boldsymbol{v}))\text{,}\]

where \(\delta(\cdot, \cdot)\) is the Kronecker delta function.

References

[Con16]

Laurent Condat. Fast projection onto the simplex and the l1 ball. Mathematical Programming, 158(1):575–585, 2016. URL: https://doi.org/10.1007/s10107-015-0946-6.

[CY14]

Jinhui Chen and Jian Yang. Robust subspace segmentation via low-rank representation. IEEE Transactions on Cybernetics, 44(8):1432–1445, 2014. URL: https://doi.org/10.1109/TCYB.2013.2286106.

[CCS10]

Jian-Feng Cai, Emmanuel J Candès, and Zuowei Shen. A singular value thresholding algorithm for matrix completion. SIAM Journal on Optimization, 20(4):1956–1982, 2010. URL: https://doi.org/10.1137/080738970.