Numerical algorithms, optimisation, and computer performance: a selection of research from 2020

en , Parallel computing , Linear algebra · 30 Dec 2020

I won’t bore you with clichés about how weird and difficult this year has been. Instead, I’d like to point out straightaway how large this year’s selection of papers is: more than 400 papers, which is more than twice the papers I selected last year.

This year, I haven’t changed my Google Scholar query:

algorithm sparse method parallel performance matrix numerical optimization graph

but I have added some more supplemental papers, as I found fit.

In the coming year, I expect to add some papers about Computational Sustainability, a field I’m getting very interested in.

Some general trends

This year, I found several papers about more or less “classic” mathematical optimisation applied to machine learning and deep learning. To me, this was quite predictable, but it’s always good to see how powerful and applicable some of those “classic” techniques are.
Multi-precision algorithms are continuing to grow in popularity in the field of numerical analysis. The idea behind these algorithms is already very in deep learning. Since we have 400 papers this year, I’m going to get some help from machine learning to analyse the data further.

An analysis of the data

Load data and fit LDA model to find the topics.

rng('default');
textData = string(fileread("bibliography.txt"));
textData = splitlines(textData);
textData = iCleanBibliographyText(textData);
documents = iPrepareDocument(textData);
[bagTraining, docsTest] = iSplitTrainingTest(documents, 0.8);

Train a number of models, changing the number of topics.

allNumTopics = 2:20;
sizeAllNumTopics = numel(allNumTopics);
mdl = cell(sizeAllNumTopics, 1);
logProbabilities = cell(sizeAllNumTopics, 1);
meanLogProbabilities = zeros(sizeAllNumTopics, 1);
perplexities = cell(sizeAllNumTopics, 1);
meanPerplexities = zeros(sizeAllNumTopics, 1);
for k = 1:sizeAllNumTopics
    rng('default');
    
    numTopics = allNumTopics(k);
    mdl{k} = fitlda(bagTraining, numTopics, 'Verbose', 0);
    logProbabilities{k} = logp(mdl{k}, bagTraining, 'NumSamples', 200);
    [~, perplexities{k}] = logp(mdl{k}, docsTest, 'NumSamples', 200);
    meanLogProbabilities(k) = mean(logProbabilities{k});
    meanPerplexities(k) = mean(perplexities{k});
end

Compare the models, plot log-probabilities over the training set and perplexities over the test set.

figure;

title("Model comparison");

subplot(1, 2, 1);

plot(allNumTopics, meanLogProbabilities);

ylabel("Log-probability on the training set");

xlabel("Number of topics");

subplot(1, 2, 2);

plot(allNumTopics, meanPerplexities);

ylabel("Perplexities on the training set");

xlabel("Number of topics");

Pick the model with the best performance on the test set.

[~, idxBestModel] = min(meanPerplexities);

bestModel = mdl{idxBestModel};

bestNumTopics = allNumTopics(idxBestModel);

iDisplayTopics(bestModel);

Topic 1: problem solution solve show order case find apply introduce approximate Topic 2: propose base many number different computational cost set good size Topic 3: numerical method first experiment accuracy fast novel two accelerate compare Topic 4: datum performance provide new require execution challenge various process feature Topic 5: model optimization learn structure network machine parameter stateoftheart present largescale Topic 6: implementation compute performance achieve gpu gpus test technique architecture design Topic 7: method problem gradient iteration complexity analysis condition optimization exist particular Topic 8: matrix memory parallel block communication distribute computation partition sparse datum Topic 9: application program code analysis performance include study hpc framework exist Topic 10: system solver large linear equation low nonlinear result example lead Topic 11: algorithm result work improve propose term vector constraint version mean Topic 12: approach time reduce work speedup well three implement experimental balance Topic 13: algorithm sparse matrix linear parallel present factorization demonstrate iterative direct Topic 14: library system paper scientific software hardware simulation highperformance spmv application Topic 15: tensor result tool paper show two decomposition factor guarantee input Topic 16: convergence function point convex rate global stochastic local bound class Topic 17: error compute precision runtime operation scale algebra computation overhead arithmetic Topic 18: graph cluster algorithm spectral random nod sequence edge quadratic match

iDisplayLogProbability(logProbabilities{idxBestModel});

Show how the best model categorises all the documents

topicMixture = transform(bestModel, docsTest);

iDisplayTopicMixtures(topicMixture, bestNumTopics);

Compare the topics found by the model

allTopicMixtures = transform(bestModel, documents);

iDisplayClosenness(allTopicMixtures);

Helper functions

function cleanDocuments = iCleanBibliographyText(document)
    % Clean the input array of strings, making sure each string
    % corresponds to a bibliography entry.
    stringSizes = strlength(document);
    
    document(stringSizes==0) = [];
    isAbstract = contains(document, "Abstract:");
    
    idxAbstract = find(isAbstract);
    % Previous abstract ends two rows before the new.
    idxEndAbstract = [idxAbstract-2; numel(document)];
    idxEndAbstract(1) = [];
    
    cleanDocuments = repmat("", [numel(idxAbstract), 1]);
    
    % Select each abstract
    for k = 1:numel(idxAbstract)
        currAbstractStart = idxAbstract(k);
        currAbstractEnd = idxEndAbstract(k);
        cleanDocuments(k) = join(document(currAbstractStart:currAbstractEnd));
    end
    
    cleanDocuments(strlength(cleanDocuments)==0 | ismissing(cleanDocuments)) = [];
end
function documents = iPrepareDocument(textData)
    % Prepare the document for analysis
    documents = tokenizedDocument(textData);
    documents = removeStopWords(documents);
    documents = erasePunctuation(documents);
    documents = removeShortWords(documents, 2);
    documents = removeWords(documents, ["Abstract", "2020", "arXiv"]);
    documents = normalizeWords(documents, 'Style', 'lemma');
end
function iDisplayTopics(mdl)
    % Display the topics in the model
    figure
    numTopics = mdl.NumTopics;
    for idxTopic = 1:numTopics
        topicWords = join(iFindTopicWords(idxTopic), " ");
        disp("Topic " + idxTopic + ": " + topicWords);
    end
    
    function words = iFindTopicWords(idxTopic)
        wordProbabilities = mdl.TopicWordProbabilities(:, idxTopic);
        wordList = mdl.Vocabulary;
        [~, idxTopProbabilities] = maxk(wordProbabilities, 10);
        words = wordList(idxTopProbabilities);
    end
end
function iDisplayTopicMixtures(topicMixtures, numTopics)
    % Display the probabily of topics by document
    figure
    area(topicMixtures);
    ylim([0 1])
    xlim([0, size(topicMixtures, 1)]);
    title("Topic mixtures over the test set")
    ylabel("Topic probability")
    xlabel("Document number")
    legend("Topic " + string(1:numTopics),'Location','northeastoutside')
end
function iDisplayLogProbability(logProbabilities)
    % Display a log-probability diagram
    figure
    histogram(logProbabilities)
    xlabel("Log probability")
    ylabel("Frequency")
    title("Document log-probabilities")
end
function iDisplayClosenness(allTopicMixtures)
    % Display the closenness between
    topicNorms = vecnorm(allTopicMixtures, 2);
    numTopics = numel(topicNorms);
    topicCloseness = (allTopicMixtures')*(allTopicMixtures)./((topicNorms')*topicNorms);
    figure    
    currCol = pcolor(topicCloseness);
    colorbar('location', 'EastOutside');
    xlabel("Topic number");
    ylabel("Topic number");
    title("Closeness between topics");
    
    plotAxis = currCol.Parent;
    plotAxis.XTick = 1:numTopics;
    plotAxis.YTick = 1:numTopics;
    plotAxis.YDir = 'reverse';
    pbaspect(plotAxis, [1, 1, 1]);
end
function [bagTraining, docsTest] = iSplitTrainingTest(documents, ratio)
    % Split training set and test set
    numDocuments = numel(documents);
    numTraining = round(ratio*numDocuments);
    documents = documents(randperm(numDocuments));
    bagTraining = bagOfWords(documents(1:numTraining));
    docsTest = documents(numTraining+1:end);
end

From the data analysis, I see:

among the list of papers, the topic of numerical methods applied to tensors occurred a couple of times (topic 15 in my analysis). This is most likely a new trend that I’d like to keep an eye on, even if I don’t know much about tensors;
reading more about the topic of error analysis for multi-precision algorithms (topic 17 in my analysis) matched what I expected from last year’s post;
the topic of matrix partitioning for parallel algorithms and distribution of matrix blocks (topic 8 in my analysis) is still occurring in my selection;
sparsity isn’t occurring that much in the papers anymore. In practice, it’s in only in one of the topics/subtopics my analysis identified (topic 13). This might be because the papers in the selection have increased, but not the papers in this area.

Abdelfattah A, Tomov S and Dongarra J (2020), "Investigating the Benefit of FP16-Enabled Mixed-Precision Solvers for Symmetric Positive Definite Matrices Using GPUs", In Lecture Notes in Computer Science. , pp. 237-250. Springer International Publishing.

[Abstract] [BibTeX] [DOI]

Abstract: Half-precision computation refers to performing floating-point operations in a 16-bit format. While half-precision has been driven largely by machine learning applications, recent algorithmic advances in numerical linear algebra have discovered beneficial use cases for half precision in accelerating the solution of linear systems of equations at higher precisions. In this paper, we present a high-performance, mixed-precision linear solver (Ax=b) for symmetric positive definite systems in double-precision using graphics processing units (GPUs). The solver is based on a mixed-precision Cholesky factorization that utilizes the high-performance tensor core units in CUDA-enabled GPUs. Since the Cholesky factors are affected by the low precision, an iterative refinement (IR) solver is required to recover the solution back to double-precision accuracy. Two different types of IR solvers are discussed on a wide range of test matrices. A preprocessing step is also developed, which scales and shifts the matrix, if necessary, in order to preserve its positive-definiteness in lower precisions. Our experiments on the V100 GPU show that performance speedups are up to 4.7× against a direct double-precision solver. However, matrix properties such as the condition number and the eigenvalue distribution can affect the convergence rate, which would consequently affect the overall performance.

BibTeX:

@incollection{Abdelfattah2020,
  author = {Ahmad Abdelfattah and Stan Tomov and Jack Dongarra},
  title = {Investigating the Benefit of FP16-Enabled Mixed-Precision Solvers for Symmetric Positive Definite Matrices Using GPUs},
  booktitle = {Lecture Notes in Computer Science},
  publisher = {Springer International Publishing},
  year = {2020},
  pages = {237--250},
  doi = {10.1007/978-3-030-50417-5_18}
}

Abdelfattah A, Tomov S and Dongarra J (2020), "Matrix multiplication on batches of small matrices in half and half-complex precisions", Journal of Parallel and Distributed Computing., 11, 2020. Vol. 145, pp. 188-201. Elsevier BV.

[Abstract] [BibTeX] [DOI]

Abstract: Machine learning and artificial intelligence (AI) applications often rely on performing many small matrix operations -- in particular general matrix-matrix multiplication (GEMM). These operations are usually performed in a reduced precision, such as the 16-bit floating-point format (i.e., half precision or FP16). The GEMM operation is also very important for dense linear algebra algorithms, and half-precision GEMM operations can be used in mixed-precision linear solvers. Therefore, high-performance batched GEMM operations in reduced precision are significantly important, not only for deep learning frameworks, but also for scientific applications that rely on batched linear algebra, such as tensor contractions and sparse direct solvers.\ This paper presents optimized batched GEMM kernels for graphics processing units (GPUs) in FP16 arithmetic. The paper addresses both real and complex half-precision computations on the GPU. The proposed design takes advantage of the Tensor Core technology that was recently introduced in CUDA-enabled GPUs. With eight tuning parameters introduced in the design, the developed kernels have a high degree of flexibility that overcomes the limitations imposed by the hardware and software (in the form of discrete configurations for the Tensor Core APIs). For real FP16 arithmetic, performance speedups are observed against cuBLAS for sizes up to 128, and range between 1.5× and 2.5×. For the complex FP16 GEMM kernel, the speedups are between 1.7× and 7× thanks to a design that uses the standard interleaved matrix layout, in contrast with the planar layout required by the vendor's solution. The paper also discusses special optimizations for extremely small matrices, where even higher performance gains are achievable.

BibTeX:

@article{Abdelfattah2020b,
  author = {Ahmad Abdelfattah and Stanimire Tomov and Jack Dongarra},
  title = {Matrix multiplication on batches of small matrices in half and half-complex precisions},
  journal = {Journal of Parallel and Distributed Computing},
  publisher = {Elsevier BV},
  year = {2020},
  volume = {145},
  pages = {188--201},
  doi = {10.1016/j.jpdc.2020.07.001}
}

Acebrón JA (2020), "A Probabilistic Linear Solver Based on a Multilevel Monte Carlo Method", Journal of Scientific Computing., 2, 2020. Vol. 82(3) Springer Science and Business Media LLC.

[Abstract] [BibTeX] [DOI]

Abstract: We describe a new Monte Carlo method based on a multilevel method for computing the action of the resolvent matrix over a vector. The method is based on the numerical evaluation of the Laplace transform of the matrix exponential, which is computed efficiently using a multilevel Monte Carlo method. Essentially, it requires generating suitable random paths which evolve through the indices of the matrix according to the probability law of a continuous-time Markov chain governed by the associated Laplacian matrix. The convergence of the proposed multilevel method has been discussed, and several numerical examples were run to test the performance of the algorithm. These examples concern the computation of some metrics of interest in the analysis of complex networks, and the numerical solution of a boundary-value problem for an elliptic partial differential equation. In addition, the algorithm was conveniently parallelized, and the scalability analyzed and compared with the results of other existing Monte Carlo method for solving linear algebra systems.

BibTeX:

@article{Acebron2020,
  author = {Juan A. Acebrón},
  title = {A Probabilistic Linear Solver Based on a Multilevel Monte Carlo Method},
  journal = {Journal of Scientific Computing},
  publisher = {Springer Science and Business Media LLC},
  year = {2020},
  volume = {82},
  number = {3},
  doi = {10.1007/s10915-020-01168-2}
}

Acer S and Aykanat C (2020), "Reordering sparse matrices into block-diagonal column-overlapped form", Journal of Parallel and Distributed Computing., 3, 2020. Elsevier BV.

[Abstract] [BibTeX] [DOI]

Abstract: Many scientific and engineering applications necessitate computing the minimum norm solution of a sparse underdetermined linear system of equations. The minimum 2-norm solution of such systems can be obtained by a recent parallel algorithm, whose numerical effectiveness and parallel scalability are validated in both shared- and distributed-memory architectures. This parallel algorithm assumes the coefficient matrix in a block-diagonal column-overlapped (BDCO) form, which is a variant of the block-diagonal form where the successive diagonal blocks may overlap along their columns. The total overlap size of the BDCO form is an important metric in the performance of the subject parallel algorithm since it determines the size of the reduced system, solution of which is a bottleneck operation in the parallel algorithm. In this work, we propose a hypergraph partitioning model for reordering sparse matrices into BDCO form with the objective of minimizing the total overlap size and the constraint of maintaining balance on the number of nonzeros of the diagonal blocks. Our model makes use of existing partitioning tools that support fixed vertices in the recursive bipartitioning paradigm. Experimental results validate the use of our model as it achieves small overlap size and balanced diagonal blocks.

BibTeX:

@article{Acer2020,
  author = {Seher Acer and Cevdet Aykanat},
  title = {Reordering sparse matrices into block-diagonal column-overlapped form},
  journal = {Journal of Parallel and Distributed Computing},
  publisher = {Elsevier BV},
  year = {2020},
  doi = {10.1016/j.jpdc.2020.03.002}
}

Aggarwal CC (2020), "The Linear Algebra of Similarity", In Linear Algebra and Optimization for Machine Learning. , pp. 379-410. Springer International Publishing.

[Abstract] [BibTeX] [DOI]

Abstract: A dot-product similarity matrix is an alternative way to represent a multidimensional data set. In other words, one can convert an n × d data matrix D into an n × n similarity matrix S = DD^T (which contains n^2 pairwise dot products between points). One can use S instead of D for machine learning algorithms. The reason is that the similarity matrix contains almost the same information about the data as the original matrix. This equivalence is the genesis of a large class of methods in machine learning, referred to as kernel methods. This chapter builds the linear algebra framework required for understanding this important class of methods in machine learning. The real utility of such methods arises when the similarity matrix is chosen differently from the use of dot products (and the data matrix is sometimes not even available).

BibTeX:

@incollection{Aggarwal2020,
  author = {Charu C. Aggarwal},
  title = {The Linear Algebra of Similarity},
  booktitle = {Linear Algebra and Optimization for Machine Learning},
  publisher = {Springer International Publishing},
  year = {2020},
  pages = {379--410},
  doi = {10.1007/978-3-030-40344-7_9}
}

Agullo E, Cools S, Fatih-Yetkin E, Giraud L, Schenkel N and Vanroose W (2020), "On soft errors in the Conjugate Gradient method: sensitivity and robust numerical detection". Thesis at: INRIA.

[Abstract] [BibTeX] [URL]

Abstract: The conjugate gradient (CG) method is the most widely used iterative scheme for the solution of large sparse systems of linear equations when the matrix is symmetric positive definite. Although more than sixty year old, it is still a serious candidate for extreme-scale computation on large computing platforms. On the technological side, the continuous shrinking of transistor geometry and the increasing complexity of these devices affect dramatically their sensitivity to natural radiation, and thus diminish their reliability. One of the most common effects produced by natural radiation is the single event upset which consists in a bit-flip in a memory cell producing unexpected results at application level. Consequently, the future computing facilities at extreme scale might be more prone to errors of any kind including bit-flip during calculation. These numerical and technological observations are the main motivations for this work, where we first investigate through extensive numerical experiments the sensitivity of CG to bit-flips in its main computationally intensive kernels, namely the matrix-vector product and the preconditioner application. We further propose numerical criteria to detect the occurrence of such soft errors; we assess their robustness through extensive numerical experiments.

BibTeX:

@techreport{Agullo2020,
  author = {Emmanuel Agullo and Siegfried Cools and Emrullah Fatih-Yetkin and Luc Giraud and Nick Schenkel and Wim Vanroose},
  title = {On soft errors in the Conjugate Gradient method: sensitivity and robust numerical detection},
  school = {INRIA},
  year = {2020},
  url = {https://hal.inria.fr/hal-02495301}
}

Agullo E, Altenbernd M, Anzt H, Bautista-Gomez L, Benacchio T, Bonaventura L, Bungartz H-J, Chatterjee S, Ciorba FM, DeBardeleben N, Drzisga D, Eibl S, Engelmann C, Gansterer WN, Giraud L, Goeddeke D, Heisig M, Jezequel F, Kohl N, Li XS, Lion R, Mehl M, Mycek P, Obersteiner M, Quintana-Orti ES, Rizzi F, Ruede U, Schulz M, Fung F, Speck R, Stals L, Teranishi K, Thibault S, Thoennes D, Wagner A and Wohlmuth B (2020), "Resiliency in Numerical Algorithm Design for Extreme Scale Simulations", October, 2020.