Thibault Clérice
email: author.thibault.clerice@chartes.psl.eu

Ground-truth Free Evaluation of HTR on Old French and Latin Medieval Literary Manuscripts

Keywords: HTR, OCR Quality Evaluation, Historical languages, Spelling Variation

As more and more projects openly release ground truth for handwritten text recognition (HTR), we expect the quality of automatic transcription to improve on unseen data. Getting models robust to scribal and material changes is a necessary step for specific data mining tasks. However, evaluation of HTR results requires ground truth to compare prediction statistically. In the context of modern languages, successful attempts to evaluate quality have been done using lexical features or n-grams.This, however, proves difficult in the context of spelling variation that both Old French and Latin have, even more so in the context of sometime heavily abbreviated manuscripts. We propose a new method based on deep learning where we attempt to categorize each line error rate into four error rate ranges (0 < 10% < 25% < 50% < 100%) using three different encoder (GRU with Attention, BiLSTM, TextCNN). To train these models, we propose a new dataset engineering approach using early stopped model, as an alternative to rule-based fake predictions. Our model largely outperforms the n-gram approach. We also provide an example application to qualitatively analyse our classifier, using classification on new prediction on a sample of 1,800 manuscripts ranging from the 9 th century to the 15 th .

Introduction

Handwritten Text Recognition (HTR) technologies have come a long way over the last five years, to the point where data mining of medieval manuscripts and HTR-supported critical editions is getting less rare nowadays, thanks in part to the user-friendliness of interfaces such as Transkribus [START_REF] Kahle | Transkribus-a service platform for transcription, recognition and retrieval of historical documents[END_REF] and eScriptorium [START_REF] Kiessling | escriptorium: an open source platform for historical document analysis[END_REF]. HTR, however, often shows limits in its ability to adapt to other scribes or periods, as it seems to fit specific scripts and languages. For example, Schoen and Saretto [START_REF] Schoen | Optical character recognition (OCR) and medieval manuscripts: Reconsidering transcriptions in the digital age[END_REF] has shown that a model trained over 1,330 lines of the 15th-century manuscript CCC 198 1 produces around 8.73% CER over test lines of the same manuscripts, drops to 14% on the same text in another manuscript from the same decade, and can go as low as 73.23% CER for a manuscript of a different text 2 even though it is at most 20 years "younger" and in the same language.

In order to evaluate the consistency of a model on an out-of-domain document such as another manuscript or a new hand, researchers usually have to create new ground-truth transcriptions to which the model predictions are compared. In this context, it seems out of reach to leverage with confidence the amount of data that remains dormant in the open data vaults of libraries such as the Bibliothèque Nationale de France (BnF) for statistical studies, making the number of 50,149 IIIF manifests catalogued by Biblissima's portal [START_REF] Frunzeanu | Biblissima's choices of tools and methodology for interoperability purposes[END_REF] promising while leaving a bitter taste of unavailability: it would require the manual transcription of at least a few hundred lines for each manuscript 3 .

To address this, we can approach this issue not as an HTR one but rather as a Natural Language Processing (NLP) task, evaluating the apparent "correctness" of the acquired text rather than its direct relationship with the digital picture of the manuscript. Evaluating new transcriptions without ground truth has been done, but mainly for OCR and non-historical documents. For modern languages, where spelling is fixed and grammar stable, a dictionary approach in combination with some n-gram statistics have provided a solid framework for establishing the probability that a document is well transcribed. However, for languages such as old French or medieval Latin, both evolving over the span of few centuries, the issue is different. For example, Camps et al. [START_REF] Camps | Corpus and models for lemmatisation and POS-tagging of old french[END_REF] has catalogued 36 forms of the word cheval (horse) in the largest available Old French corpus. A Dictionary approach would already prove to be complex, but to make things worse, the abbreviated nature of medieval texts would require taking into account several abbreviation systems, making it unsustainable.

HTR is most often, in the humanities, not a task in itself but rather a preliminary step for corpus building (such as digital editions) or corpus analysis. In this context, HTR quality can be of primordial importance, depending on the task at hand. While Eder [START_REF] Eder | Mind your corpus: systematic errors in authorship attribution[END_REF] has suggested that good classification in stylometry is still possible for corpora with noise levels as high as 20%, even for the smallest feature sets, Camps et al. [START_REF] Camps | Handling heavily abbreviated manuscripts: HTR engines vs text normalisation approaches[END_REF] demonstrated that, for HTR, noise leads to accumulating errors throughout its post-processing (word segmentation, abbreviation resolution, lemmatization and POS-tagging), making the post-processed textual features less reliable than original character n-grams. For some other tasks, such as in corpus linguistics (e.g. semantic drift studies), the study of abbreviation systems such as the one performed by Honkapohja and Suomela [START_REF] Honkapohja | Lexical and function words or language and text type? Abbreviation consistency in an aligned corpus of Latin and Middle English plague tracts[END_REF] or even the training of large language models such as MacBerth [START_REF] Manjavacas | Adapting vs. Pre-training Language Models for Historical Languages[END_REF] might require a higher level of precision.

As such, evaluating the textual quality of an automatic transcription "from afar" is extremely useful, as it provides solid grounds to either exclude documents from analysis or help guide ground-truth creation campaigns in well-funded projects. For cultural heritage institutions, it can also provide a welcome indicator for the document that could be ingested by a research engine. We can even imagine situations where these institutions transcribe only a sample of each element of their collection, and only fully and automatically transcribe the ones that reach a certain level of quality, thus saving energy and ultimately budget on the computation front.

From a human reader's perspective, Springmann et al. [START_REF] Springmann | Automatic quality evaluation and (semi-) automatic improvement of OCR models for historical printings[END_REF] and Holley [START_REF] Holley | How good can it get? analysing and improving ocr accuracy in large scale historic newspaper digitisation programs[END_REF] have set a limit of a CER below 10% for a good OCR quality. Recently, Cuper [START_REF] Cuper | Examining a multi layered approach for classification of ocr quality without ground truth[END_REF] has proposed the evaluation of OCR quality for heritage text collections, specifically Dutch newspapers from the 17 th century, to distinguish good OCR from bad, using the aforementioned threshold. They provide a tool, QuPipe, which offers binary classification capacities, putting text in either the range [0; 10]% of CER or in the remaining range of "bad" OCR. In 2022 as well, Ströbel et al. [START_REF] Ströbel | Evaluation of htr models without ground truth material[END_REF] addressed this issue regarding HTR of cultural heritage documents, specifically from the 16 th century. They provide a strong argument for using lexical features and (pseudo-)perplexity scores for HTR quality estimation, with the specific limitation that the texts they studied, 16th-century Latin correspondence, does not provide as much variation as older languages such as historical German. We also note that correspondence may be less abbreviated, and that this dataset spans a very short period.

In this paper, we address this issue as a supervised classification task, based on a dataset of around 50,000 lines of ground truth spanning from the 9th through to the 15th century. Following the conclusion of Cuper [START_REF] Cuper | Examining a multi layered approach for classification of ocr quality without ground truth[END_REF], we augment the number of categories we want to find: we distinguish Good ([0, 10)%), Acceptable ([10, 25)%), Bad ([25, 50)%), and Very Bad (≥ 50%) rates of OCR. This provides a more fine-grained evaluation of the transcription and allows for guided transcription campaigns, by addressing either the low-hanging fruits (Acceptable) or the rotten ones. We evaluate three kinds of basic architectures (GRU with attention, BiLSTM and TextCNN) on line classification using real-life "bad" transcriptions and precomputed CER scores.

The resulting models have shown promising results, with quality levels such as Very Bad and Good being well recognized. In order to evaluate the models and showcase their usefulness, we also provide an example of a real-life classification application, where 1800 manuscripts were randomly selected from the BnF and classified by our best model.

In summary, the contributions of this paper are:

1. a new approach for HTR evaluation of historical languages with variable spellings; 2. a new method to produce ground truth for OCR evaluation that does not rely on artificially and manually tuned generation; 3. an initial evaluation of the output and a quick glance at the state of HTR for Old French and Medieval Latin over six centuries.

The remainder of this paper is organised as follows. We start by addressing the background in Section 2, specifically regarding the specifics of Old French and Medieval Latin and the idea of readability. In Section 3, we describe the HTR datasets we used and their particularities. In Section 4, we describe the architecture of the models, their feature engineering and the process behind the generation of bad predictions. In Section 5, we describe the set-up of our model selection and evaluation. Finally, in Section 6, we analyse the result both on the dataset produced ad hoc (described in Sections 3 and 4), but also on completely unseen documents from the BnF, to showcase the capacities of such models.

Background and Related Work

Handwritten Text Recognition, a sibling or sub-task of Optical Character Recognition, aims at recognising text from digitised manuscripts. In the last five years, the digital humanities landscape has seen a surge in HTR engines, as well as transcription interfaces that connect and work well with these engines, from the dominant Transkribus [START_REF] Kahle | Transkribus-a service platform for transcription, recognition and retrieval of historical documents[END_REF] to the open-source pair of eScriptorium [START_REF] Kiessling | escriptorium: an open source platform for historical document analysis[END_REF] and Kraken [START_REF] Kiessling | The Kraken OCR system[END_REF]. To be able to recognize text, users have to provide models, which are themselves the result of supervised training on ground truth data (human provided transcriptions).

Printed books have been, over the last few decades, the focus in terms of remediation, from their analogue form to a digitized picture and finally to a machine-readable (and human searchable) text. With the advances in HTR over the last five years, the focus can now shift or be shared with materials that have, for the most part, remained inaccessible from a digital point of view, except for pictures. Latin manuscripts are present during the whole period of manuscript production in western Europe. Literary Old French manuscripts exist from the 12 th century onward, with only a hundred known surviving manuscripts in the 12 th century [START_REF] Careri | Livres et écritures en français et en occitan au XIIe siècle: catalogue illustré[END_REF]. Over the span of these seven centuries, multiple forms of handwritten scripts have existed, for both French and Latin. As an example, the 2016 ICHFR Competition on the Classification of Medieval Handwritings in Latin Script [START_REF] Cloppet | ICFHR2016 competition on the classification of medieval handwritings in latin script[END_REF] provided ground truth for the classification of 12 main families, of which at least six are represented in our datasets. This diversity makes training models for HTR quite complex but also a reachable goal, as they tend, specifically for literary manuscripts, to be more readable and stable between different hands.

Medieval French and Latin present both dialectal and scriptural variation in synchrony on top of diachronic evolution. Old French's syntax varies chronologically and geographically. The spelling is simply variable. While Latin shows some level of variation, it differs from Old French mostly in its higher rate of abbreviation. These observations are limited to the context of the datasets at hand, which are literary works (including scholastic, theological and medical works). The Old French CREMMA Medieval dataset [START_REF] Pinche | Cremma medieval[END_REF] has 0.97% of horizontal tildes and 0.16% of vertical ones, which are markers used in the dataset guidelines to indicate various similar abbreviation diacritics [START_REF] Pinche | Guide de transcription pour les manuscrits du Xe au XVe siècle[END_REF]. Using the same guidelines, the CREMMA Medieval Latin dataset shows a rate of 5.63% and 1.52% for the same characters. This difference could be due to the nature of the transcribed texts.

The question of abbreviation and the specificity of medieval literary manuscripts has provoked many discussions in terms of how to transcribe documents, from a completely "diplomatic" approach with variants of letters to "semi-diplomatic" approaches. In the last year, three authors have provided guidance or thoughts around guidelines for transcriptions: Pinche [START_REF] Pinche | Guide de transcription pour les manuscrits du Xe au XVe siècle[END_REF] focusing on Old French, Schoen and Saretto [START_REF] Schoen | Optical character recognition (OCR) and medieval manuscripts: Reconsidering transcriptions in the digital age[END_REF] on Middle English, and Guéville and Wrisley [START_REF] Guéville | Transcribing medieval manuscripts for machine learning[END_REF] on Latin. The CREMMA guidelines have been used by 5 other datasets for a total of 1.15 millions of characters over fifty manuscripts, which make them the most diverse and comprehensive ones for HTR of medieval manuscripts in Latin and Old French.

The most traditional metrics for HTR and OCR are both Word Error Rate (WER) and Character Error Rate (CER). The first one proves to be complicated to apply in Old French and Medieval Latin, as spaces in medieval manuscripts tend to vary in size or simply be nonexistent from a modern perspective, relying on the knowledge of the reader to separate words-or the ability of NLP models to separate them [START_REF] Clérice | Evaluating Deep Learning Methods for Word Segmentation of Scripta Continua Texts in Old French and Latin[END_REF]. The second one works well, with the limitation that spaces are often the first source of mistakes. CER corresponds to the sum of character insertion, removal and replacement over the total number of characters, thus providing a finegrained metric.

As mentioned earlier, in the introduction, both CER and WER require ground truth, and other metrics currently discussed as alternatives, such as the (pseudo-)perplexity or lexical measures proposed by Ströbel et al. [START_REF] Ströbel | Evaluation of htr models without ground truth material[END_REF]. The other approach to evaluating quality without ground truth is to predict a class of CER, such as the work done by Bazzo et al. [START_REF] Bazzo | Assessing the impact of OCR errors in information retrieval[END_REF]. These approaches rely on features such as n-grams, word statistics and language classifier outputs which are difficult to leverage in the present context. In order to train their classifier, Bazzo et al. [START_REF] Bazzo | Assessing the impact of OCR errors in information retrieval[END_REF] and Nguyen et al. [START_REF] Nguyen | Readocr: A novel dataset and readability assessment of ocred texts[END_REF] engineered bad predictions by creating rules to reproduce the most common errors in OCR, such as "rn" becoming "m". These bad predictions are then fed to their model along with the metrics both papers want to predict.

Nguyen et al. [START_REF] Nguyen | Readocr: A novel dataset and readability assessment of ocred texts[END_REF] provide an innovative approach to the issue of noise in OCR by shifting from a CER/WER problem to a readability one: if the reader "can reod a txt with miffpelling" without having to refer back to the picture, at least one of the goals of OCR has been achieved. As simply put by Martinc et al. [START_REF] Martinc | Supervised and unsupervised neural approaches to text readability[END_REF], "Readability is concerned with the relation between a given text and the cognitive load of a reader to comprehend it". It is even more important in the context of handwritten documents where a somewhat badly but readable HTR output can be easier for non-specialists to read than the original. In the field of readability assessment, Martinc et al. [START_REF] Martinc | Supervised and unsupervised neural approaches to text readability[END_REF] has shown that supervised models perform adequately, while Nguyen et al. [START_REF] Nguyen | Readocr: A novel dataset and readability assessment of ocred texts[END_REF] has shown that this translates to the OCR issues as well. This has not been applied to any medieval dataset that we know of.

Dataset

To train different models, we reused the data from various projects, aligned with the same guidelines used by Pinche [START_REF] Pinche | Guide de transcription pour les manuscrits du Xe au XVe siècle[END_REF]. Our experiment was made possible by the open release of many projects' datasets, including one MA thesis and one student project [START_REF] Vlachou-Efstathiou | lat.o.41 -eutyches "de uerbo[END_REF][START_REF] Biay | Tnah-2021-decameronfr[END_REF]. We used the ground truth of the CREMMA [START_REF] Clérice | Cremma manuscrits médiévaux latins[END_REF] and CREMMALab [START_REF] Pinche | Cremma medieval[END_REF] projects, the Rescribe [START_REF] White | Caroline minuscule by rescribe[END_REF] project, and the GalliCorpora [START_REF] Pinche | Données HTR incunables du 15e siècle[END_REF][START_REF] Pinche | Données HTR manuscrits du 15e siècle[END_REF] projects, for a total of 42,292 lines (see Table 1). We include one dataset of incunabula, which use graphical shapes similarly to literary manuscripts (but with more regularity), while also using an abbreviation system.

Table 1

Training material for our models and our future bad transcription dataset.

The datasets present not only two main languages but also many different levels of digitization quality (including old binarization), different kinds of handwriting families, different abbreviation levels and different genres. For example, while the CREMMA Medieval dataset focuses more on literary texts, specifically hagiographical and chanson de geste texts, the CREMMA Medieval LAT corpus offers theological commentaries and medicinal recipes, each genre having its own specific vocabulary. The dataset in general is skewed towards French and the gothica handwritten family.

The transcription guidelines of Pinche [START_REF] Pinche | Guide de transcription pour les manuscrits du Xe au XVe siècle[END_REF] provide simplification rules: allographic approaches are forbidden (different shapes of s such as long s and "modern" s are not differentiated), macrons and general horizontal-line diacritics over the letters such as tildes are represented by horizontal tildes, any "zigzag" 4 or similarly shaped forms are simplified into superscript vertical tildes, etc. This allows for simpler transcriptions and also limited diversity of characters for the machine to learn, satisfying both the human transcriber in terms of the learning curve of the guidelines, and the HTR engine in terms of complexity. Each corpus was passed through the ChocoMufin software [START_REF] Clérice | Choco-Mufin, a tool for controlling characters used in OCR and HTR projects[END_REF] using project-specific character translation tables. This software, along with these tables, allows each dataset to be controlled at the character level and adapted to guideline modifications. It also allows project-specific transcription standards to be translated to a more common one, such as Pinche's.

Proposed Method

Our goal is to be able to predict a quality class for any HTR output on medieval French and Latin. First, we design a way to generate ground truth for the quality assessment of HTR output. Then, we propose three supervised text-based models, with specific adaptations to handle both languages with a single classifier.

"Bad Prediction" Ground Truth

In order to train our classification model, we require ground truth material along a CER class: Good ([0; 10)%), Acceptable ([10; 25)%), Bad ([25; 50)%) and Very Bad (≥ 50%). In order to have real life errors, and to reproduce the rather difficult to predict capacity of a model to confuse certain characters with others in specific settings, we propose a three-step method:

1. We train Kraken [START_REF] Kiessling | The Kraken OCR system[END_REF] models based on the complete dataset, or on a subset. We voluntarily stop some of the training in very early stages, when the CER on the validation dataset remains high. We also keep one "best" model [START_REF] Clérice | Generic CREMMA Model for Medieval Manuscripts (Latin and Old French)[END_REF] trained on the full dataset. 2. We run each model on our two biggest and most diverse repositories, CREMMA Medieval and CREMMA Medieval LAT. We also run a model trained on modern and contemporary scripts, Manu McFrench [START_REF] Chagué | HTR-United -Manu McFrench V1 (Manuscripts of Modern and Contemporaneous French)[END_REF] to create garbage-level transcriptions. 3. We evaluate each line's CER and store it alongside the line. We also keep the ground truth, whose CER is estimated as 0. We remove short lines (fewer than 15 characters) and duplicated predictions across models for the same line.

Regarding the final models for prediction production, we have 16 models, allowing for a maximum of 16 versions of each line, if none of the models predict the same text (see Table 2 for examples):

1. 4 models trained on the same train and validation dataset as best with a validation CER of 55.9, 28.3, 23 and 20.8% according to Kraken. 2. 5 models trained on the CREMMA Medieval LAT dataset only, from the 1 st to the 6 th epochs, ranging from 86% to 46% of CER. These provide variable CER on unseen data from the test set of both CREMMA dataset but also on training and validation sets as they did not reach their full capacities during the training phase. After filtering small and repeated predictions, we have access to 322,903 lines of "HTR Predictions, CER" couples (see in appendix Figure 6). We then translate that into each bin of CER to produce the four established classes.

Model Architecture

We applied three model architectures, common to many NLP task, with an embedding-sentence encoder-linear classifier structure where only the sentence encoder changes from one model to another (see Figure 2). The embedding layer takes into account special tokens (Padding, Unknown char, Start of Line, End of Line) and each character according to the Unicode NFD5 normalization of the line, for which characters and their diacritics are separated, e.g. becomes

[e]+[´].
The linear layer is a simple (Encoding Output Dimension, Class Count) decision layer. Each model uses a cross-entropy loss function 6 and reduces its learning rate at plateau using the validation set's macro averaged recall metric. Optimization of the model is done through the Ranger optimizer [START_REF] Wright | New deep learning optimizer, ranger: Synergistic combination of RAdam + LookAhead for the best of…[END_REF].

The encoding layer varies between three different forms:

• The first version uses a single BiLSTM network where the sentence encoding is the result of the concatenation of the start-of-line token (BOS) and end-of-line token (EOS) hidden state. • The second version follows the architecture of sentence-level attention proposed by Yang et al. [START_REF] Yang | Hierarchical attention networks for document classification[END_REF], using a bidirectional GRU. The encoded sentence vector is the sum of products of the hidden state of each token with its attention. Attention is also provided as an output for human interpretation of the results. • The last one, TextCNN [START_REF] Kim | Convolutional neural networks for sentence classification[END_REF], uses the concatenation of the Max Pooling of each n-gram size (2, 3, 4, 5, 6) taken into account by a convolutional neural network.

As we deal with two different languages, we added another special token, following the work of Martin et al. [START_REF] Martin | Controllable Sentence Simplification[END_REF] and Gong et al. [START_REF] Gong | Enriching word embeddings with temporal and spatial information[END_REF]: for each encoding variation we add one variation of the codec where the first token after the beginning-of-string is a metadata token indicating the language. Thus, a line such as Fra on de q ̃l vertuses petis pies sont que vo will be encoded as <BOS><FRO>Fra on de q ̃l vertuses petis pies sont que vo<EOS>.

Experimental Setup

In order to avoid lexical bias and to ensure the strength of our analysis, we propose a 5-Foldlike experiment, where each subset for train, validation and test are the results of split across manuscripts. For each K, two French manuscripts and two Latin ones are used for the validation set and the test set, and they differ by at least one manuscript from one K to another, leaving three K completely different (K1, K3, K5; see Table 3). Each test set also contains a Latin manuscript that was not used in any of the HTR model training or validation: Berlin, Hdschr. 25. This manuscript was then used for model evaluation, to have a stable pillar for evaluation. Models are then evaluated using class-specific precision and recall, as well as macro averaged precision and recall.

For our baseline, we use the relative frequency of the 2000 most common n-grams of size 3, 4 and 5 as features and feed them to a linear classifier, with cross entropy loss and the Adam optimizer. We run each model architecture once for each K, resulting in 7 different results with the baseline (presence/absence of language token for the three encoding modules + baseline).

Our whole pipeline uses pandas for data preparation [START_REF] Mckinney | pandas: a foundational python library for data analysis and statistics[END_REF], PyTorch [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF] for model development, and Pytorch Lightning [START_REF] Falcon | The PyTorch Lightning team, PyTorch Lightning[END_REF] for the training, evaluation and prediction wrapping.

Table 4

Test results statistics for each K and each model configuration.

Experiments

Model Classification Results

The first conclusion we can draw from the experience is that our models always beat the baseline (see Table 4 and, in the Appendix, Table 4 for more details). No RNN-based architecture clearly beats the other, but TextCNN clearly underperforms. The introduction of the language metadata token helps when detecting Good transcriptions (delta ≈ +7% for attention's median precision, ≤ +1% for the recall) for both RNN based models. Models without a language marker tend to outperform models with language markers, except for the Very Bad class where the delta is up to +6% in favour of models without language tokens (using median precision scores).

Regarding the variability of results, we found that the length of the string had an impact on the prediction, no matter the model architecture. Surprisingly, none of the models withstand long noisy lines: the accuracy of the Very Bad class is inversely correlated with line size. On the contrary, depending on the encoder, some classes benefit from longer strings: Good lines benefit from it with all models except the baseline. TextCNN is the only model to really correlate accuracy on the Bad and Acceptable classes with line length.

Finally, for all models except the baseline, the most common confusion is always in the "adjacent" class(es) (see Figure 4). For the classes Acceptable and Bad, which have two neighbours, the error rate is evenly split between them: the class Acceptable tends to be confused with either Good or Bad. This shows the model's ability to understand cleanness or noise, but also shows the limit of these classes: for a line with 50 characters, such as "quãt tel eufaut gist en tes lieu. Derite respoint", 6 mistakes are enough to swing into the Acceptable category (Ground truth: "quãt tel enfant gist en tel lieu . Uerite respon", one space has been removed before the dot).

Overall, with an accuracy for the Good and Very Bad classes around 50% on these languages, and considering that most of the confusions are from adjacent classes (e.g. Good is confused with Acceptable, Acceptable with Good and Bad, etc.), the solution performs well either at filtering badly read manuscripts, or keeping only the very good ones. The Acceptable class and the Bad class have stable performance facing variable line length, although the Acceptable class shows the worst classification performance.

Application on a Real-World Library Dataset

As a real-world application, we wanted to apply one of our best models to an unseen dataset, in the same way that we envision cultural institutions might use the tool. We describe the set-up for this particular experiment below, and then evaluate the results of the classification model with regard to the capacity of the HTR model; we also study some randomly sampled elements.

Set-up

To evaluate on as much unseen data as possible, we crawled the Biblissima IIIF collection portal [START_REF] Frunzeanu | Biblissima's choices of tools and methodology for interoperability purposes[END_REF]. We searched individually for each combination of language (French, Latin) and century (9 th to 15 th ), limiting the number of samples retrieved to 500 manuscripts. We then sam- pled 10 sequential pictures from each manuscript. 7 To avoid empty pages (which tend to be at the start and the back of each book's digitization or IIIF manifest at the BnF), we take either the ten first pictures from the second decile of the manifest, or from the 20 th up to the 30 th if there are fewer than 100 pictures, or the 10 last if there are fewer than 20 pictures.

Each downloaded sample is then segmented using YALTAi [START_REF] Clérice | You Actually Look Twice At it (YALTAi): using an object detection approach instead of region segmentation within the Kraken engine[END_REF] with the included model designed for cultural heritage manuscripts and the base Kraken BLLA segmenter [START_REF] Kiessling | A modular region and text line layout analysis system[END_REF]. As YAL-TAi provides different zones-from the margin to main body of text-through numbering, we only consider lines that are part of the main bodies of text of each model, thus excluding any marginal or paratext. We then use Kraken to predict a transcription for each line, using the best trained model as described in our first experiment. Next, we feed each line to our best BiL-STM model (K-Fold 1 has the best recall/precision on Good) while keeping the line metadata: language, century, manuscript identifier, and page identifier.

Finally, we provide three different evaluations of the transcriptions. The first is based strictly on the number of lines predicted in each class (Good, Acceptable, etc.). The second is page-based: we take the most common prediction for all lines. The last one is manuscript based: we take the most common page prediction, using the previous page-based metric.

Evaluation

Overall, the HTR prediction results produced by our BiLSTM module are in line with the HTR strength on the dataset (see Figure 5). The model performs extremely well on early manuscripts thanks to the presence of two datasets of early manuscripts (Eutyches and Caroline Minuscule) It performs well on Old French except for the 13th century, where Bad predictions are more common. The relative frequency of Very Bad predictions tends to grow as we get closer to the 16th century: from the data we have seen, this could be due to the presence of non-literary manuscripts written in cursive, for which our model has no ground truth.

If we look at the sampled predictions (Appendix, Table 2), most Good predictions seem correct or nearly correct. However, we can see that the metadata from Biblissima and the BnF has some limitations when used automatically, as it can produce problematic results: most 12 th century Acceptable predictions are probably in Latin, which would indicate a multilingual manuscript or a badly catalogued one. This issue also arises in the crawler for the century, as some manuscripts were catalogued as French but with a production date that is before the first known French document: these are most likely multilingual documents, with either a collection of various leaves from previous manuscripts, or the inclusion of the language used for marginal notes. 3 out of the 6 Acceptable predictions between the 13 th and the 14 th century are definitely readable and understandable, and we cannot but wonder if the lack of spaces in "q̃merueilles fu lacitebiengarne mlt" is responsible for its classification as Acceptable rather than Good. We note that at least one Very Bad prediction in French, "OU EtE L. Cheualier de Monifort, son Oncle, Gles", seems rather readable, albeit with more corrections than for a Good transcription. Latin shows the same trend, in being accurate over Good and Acceptable.

Conclusion

The ability to filter, without pre-transcribing samples, automated transcriptions of manuscripts in Latin, Old French or any other Western historical language, might lead to the production of datasets designed for analysis that relies on better transcriptions, or to guiding cultural heritage institutions and their partners in the production of new ground truth. Producing HTR ground truth does indeed require time, skilled transcribers and, last but not least, budget. However, most current error rate prediction or HTR output analysis models rely on n-gram frequencies and lexical features-two approaches that are often less viable for languages such as Old French which "suffers" from a highly variable spelling system or for languages like Latin which are potentially highly abbreviated, with abbreviations changing even within a single manuscript, depending on the context, the topic and the scribe.

In this context, we chose to treat CER range predictions as a sentence-like classification problem, for which we implemented three basic models, using either a single BiLSTM encoder, an attention-supported GRU, or a TextCNN encoder. These three tools show stronger results than an n-gram based baseline. On top of this, we include a language metadata token which can improve the reliability of the lowest range of CER (between 0 and 10%, the Good class) while worsening the classification's reliability for the highest range (over 50%, the Bad class). For the purpose of training these models, we propose a new way to generate real life "bad transcriptions", using early-stopped HTR models, or models trained on small samples of data: this provides an alternative to previous rule-based generation of "bad transcription" ground truths.

We show that on a completely unknown dataset of around 1,800 manuscripts, analysed with a new HTR model specifically trained on medieval Latin and French, the number of welltranscribed manuscripts predicted is on par with the ground truth for that dataset. The quality assessment predictions provide quick insights for larger collections, and could be run relatively often by cultural heritage institutions.

In the future, hyper-parameter fine-tuning and other encoders could be used in the architecture. Specifically, with more correctly transcribed manuscripts, including the abbreviations in their transcriptions, fine-tuning larger language models could allow the application of (pseudo-)perplexity ranking such as the one proposed by Ströbel et al. [START_REF] Ströbel | Evaluation of htr models without ground truth material[END_REF], while allowing for partial noise in the training data. We hope to see such classification of manuscripts used by ground truth producers in order to enhance the robustness of openly available HTR models.

A. Appendix

The software has been archived at the following address: https://doi.org/10.5281/zenodo. 7233984. A good chunk of the data is available here: https://github.com/PonteIneptique/ neNequitia/releases/tag/chr2022-release.

Manuscripts metadata and the predictions in XML ALTO formats for Section 6 are available at https://doi.org/10.5281/zenodo.7234399. The same repository contains also the XML data for training the classifier.

Figure 1 :

 1 Figure 1: Example of lines. (a) comes from the GalliCorpora manuscript dataset, (b) from the incunabula one. (c) is drawn from the Eutyches MA Thesis, (d) and (e) from the CREMMA Medieval (French) dataset. (f) and (g) are both taken from the CREMMA Medieval Latin repository.

  3. 1 model trained on the Eutyches (Latin, Carolingian of the 9 th century) and the Decameron (French, 16 th century) datasets with a 98.5% CER on its validation set. 4. 3 models trained on the CREMMA Medieval (Old French) dataset only, fine-tuned from the Manu Mc French Model, from 11% of CER down to 8.2%. 5. Manu McFrench, the best model and the ground-truth data.

Figure 2 :

 2 Figure 2: Available model architectures. Elements in orange are optional or varying elements, elements in blue are common to all models.

K

Figure 3 :

 3 Figure 3: Regression of accuracy based on lines' length overall 5-Fold test sets. Common manuscript not included (Berlin, Hdschr. 25).

Figure 4 :

 4 Figure 4: Confusion rate dispersion in the errors made by each model. Only confusion that happens more than 50 times is taken into account, as well as the total number of errors greater or equal to 300. The graph can be read as follows: for the baseline, 40% of the errors for the ground truth class Good are Acceptable predictions.

Figure 5 :

 5 Figure 5: Predictions distribution per line (first two rows), per page (row 3 & 4), per manuscript (last rows) over languages and centuries filtering.

Figure 6 :

 6 Figure 6: "Bad transcriptions" CER Violin plot, per manuscript. Most manuscript have a strong enough diversity of CER to train upon.

Table 2

 2 Example of pairs of predictions for the same line for a file of CREMMA Medieval (University of Pennsylvania 660, Le Pélerinage de Mademoiselle Sapience). The first line is the ground truth, the second our best model trained on the full dataset for production, the 4 th from the bottom is from Manu McFrench. Note that the diacritics are not consistently transcribed.

	[é]

  BnF fr. 25516 BnF fr. 3516, BnF fr. 25516 BnF fr. 24428, BnF. Arsenal 3516 BnF fr. 24428, BnF fr.844 Pennsylvania Codex 909, BnF fr.844

			1	2	3	4	5
	Validation French BnF fr. 17229, Latin Arras 861, CCCC Ms 165	CLM 13027, CCCC 165	CLM 13027, Montpellier, H318	BnF lat. 6395, Montpellier, H318	BnF lat. 6395, Laur. Plut.33.31
	Test BnF Arsenal 3516Train French Good	80,056	76,564	65,764	39,165	39,165
		Acceptable	44,346	41,769	34,429	35,803	35,803
		Bad	60,381	59,265	51,637	41,793	41,793
		Very Bad	71,008	71,053	60,898	52,212	52,212
	Validation Good	4,246	98,57	12,770	11,625	11,625
		Acceptable	3,933	10,377	12,496	8,492	8,492
		Bad	4,338	10,884	13,430	10,250	10,250
		Very Bad	4,867	15,428	18,386	11,461	11,461
	Test	Good	9,165	7,046	14,933	42,677	42,677
		Acceptable	9,744	5,877	11,098	13,728	13,278
		Bad	12,763	7,333	12,415	25,439	25,439
		Very Bad	18,056	7,350	14,647	30,258	30,258

, BnF fr. 13496 BnF fr. 24428, BnF fr. 411 BnF fr. 844, BnF fr. 22549 BnF fr.412, Phil., Col. of Phys. 10a 13 Bodmer 168, Vat. reg. lat. 1616 Latin Sorbonne Fr. 193, CLM 13027 CCCC Ms. 236, H318 BnF lat. 6395, Egerton 821 BnF fr. 16195, Laur. Plut. 33.31 Laur. Plut. 53.08, BnF lat. 8236

Table 3

 3 Composition of K-Folds set, based on manuscript selection.

Table 6 :

 6 Examples of HTR Prediction on unseen documents and their classification by the model.

	Lang Century Prediction	Transcription
	fro	12 Good	ra monstre de couf de uoudenay
	fro	12 Good	uucuns pair er que Iehuz de le chaulre le ieuue uor de teuteur dicelui
			office oy a este priuoz et de loucez pour msen et acaise de cereus cu ⁊
			decpa
	fro	12 Good	seriant estoit exilliez en laueniance de sacolpe. li poures
	fro	13 Good	les cõdurroit car il sauoit crop bien. coz lespas. ⁊
	fro	13 Good	Se il lẽ set dire nouele
	fro	13 Good	tiseras. ⁊ ieres sempres amendes. ⁊ en un au
	fro	14 Good	Procureur du Roi du mẽme iour qui ne l empeche. lOrdon.
	fro	14 Good	sonpere tous les rodais et les tartcites

Five million lines would be required for the mentioned set of manifests of the BnF with only 100 lines per manuscript. As a comparison point, the accumulated number of lines of manuscript dataset, regardless of the script or language, publicly available on the HTR-United catalog[START_REF] Chagué | HTR-United: Ground Truth Resources for the HTR and OCR of patrimonial documents[END_REF] is 164,418 at the end of August 2022.

Official name from the Unicode specifications for the character U+299A.

Normalization Form Canonical Decomposition.

Code available at https://github.com/PonteIneptique/neNequitia.

Note that we are not talking about pages but about pictures: in some cases, most commonly in the case of digitised microfilms, one picture can contain two pages.

Acknowledgments

I want to thanks Jean-Baptiste Camps, Ariane Pinche and Malamatenia Vlachou-Efstathiou for their constant feedback and replies on some particular questions regarding manuscripts or HTR data. Thank you to Ben Nagy for his proof-reading of the pre-print version.

This work was funded by the Centre Jean Mabillon and the DIM MAP (https://www. dim-map.fr/projets-soutenus/cremmalab/).


Oleg Balabanov

Matthias Beaupère

Laura Grigori

Victor Lederer

Block subsampled randomized Hadamard transform for low-rank approximation on distributed architectures

Keywords: randomization, sketching, embedding, low-rank approximation, parallel computing

This article introduces a novel structured random matrix composed blockwise from subsampled randomized Hadamard transforms (SRHTs). The block SRHT is expected to outperform well-known dimension reduction maps, including SRHT and Gaussian matrices, on distributed architectures with not too many cores compared to the dimension. We prove that a block SRHT with enough rows is an oblivious subspace embedding, i.e., an approximate isometry for an arbitrary low-dimensional subspace with high probability. Our estimate of the required number of rows is similar to that of the standard SRHT. This suggests that the two transforms should provide the same accuracy of approximation in the algorithms. The block SRHT can be readily incorporated into randomized methods, for instance to compute a low-rank approximation of a large-scale matrix. For completeness, we revisit some common randomized approaches for this problem such as Randomized Singular Value Decomposition and Nyström approximation, with a discussion of their accuracy and implementation on distributed architectures.

Introduction

Randomization has become a powerful tool for tackling massive problems in numerical algebra and data science [START_REF] Martinsson | Randomized numerical linear algebra: Foundations and algorithms[END_REF][START_REF] David | Sketching as a tool for numerical linear algebra[END_REF][START_REF] Vershynin | High-dimensional probability: An introduction with applications in data science[END_REF][START_REF] Michael | Randomized algorithms for matrices and data[END_REF]. Modern randomized methods can, in particular, provide solutions to problems of dimensions beyond the reach of deterministic methods, and allow effective use of computational resources. Recent significant development has made them very reliable, and not just used as a last resort, as it was not so long ago. Along with increased efficiency, they can now provide strong accuracy guarantees with a user-specified failure probability that can be chosen extremely low, say 10 -10 , without much impact on computational costs.

This article is concerned with randomized methods that are based on a dimension reduction, called sketching [START_REF] David | Sketching as a tool for numerical linear algebra[END_REF], with oblivious ℓ 2 -subspace embeddings (OSEs) defined below. Definition 1.1. Let 0 ≤ ε < 1 and 0 < δ < 1. A random matrix Ω ∈ R l×n is said to be a (ε, δ, d) OSE, if for any fixed d-dimensional subspace

V ⊆ R n , ∀x ∈ V, |∥x∥ 2 2 -∥Ωx∥ 2 2 | ≤ ε∥x∥ 2 2 (1.1)
holds with probability at least 1 -δ.

It is a consequence of the Johnson-Lindenstrauss lemma [START_REF] William | Extensions of lipschitz mappings into a hilbert space 26[END_REF] that there exist (ε, δ, d) OSEs of sizes l = O(ε -2 (d + log 1 δ ). The fact that n does not appear in the right-hand-side and the logarithmic dependence on the

Contributions

The proposed block SRHT matrix has the potential to combine the benefits of structured oblivious embeddings, such as the SRHT, with the benefits of unstructured ones, such as Gaussian, from the complexity and performance standpoint. We are not aware of any other work in this direction. We prove that the block SRHT matrix of size l = O(ε -2 (d + log n δ ) log d δ ) satisfies the (ε, δ, d) OSEs property. This result is similar to that for standard SRHT from the literature [START_REF] Tropp | Improved analysis of the subsampled randomized Hadamard transform[END_REF][START_REF] Boutsidis | Improved matrix algorithms via the subsampled randomized hadamard transform[END_REF] and implies that the two matrices should yield approximations of similar quality. This result does not simply follow from the analysis in [START_REF] Tropp | Improved analysis of the subsampled randomized Hadamard transform[END_REF][START_REF] Boutsidis | Improved matrix algorithms via the subsampled randomized hadamard transform[END_REF], and particularly requires incorporation of a useful technical trick that, to the best of our knowledge, was not employed before in the randomized numerical linear algebra community.

In addition, we present low-rank approximation methods such as RSVD [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF], Nyström approximation, and the single-view approximation from [START_REF] Tropp | Streaming low-rank matrix approximation with an application to scientific simulation[END_REF] in a new unified projection-based form that clearly shows their connection. We provide a rigorous characterization of their accuracy based solely on the OSE property, and therefore compatible with any types of sketching matrices that satisfy this property, including the block SRHT. For RSVD and Nyström approximation, this characterization follows almost directly from standard results from the literature. On the other hand, for the single-view approximation, our results are new. In particular, the novelty lies in the use of the projection-based interpretation of the single-view approximation to show that it is almost as accurate as RSVD if the embeddings for the "core sketch" are OSEs of sufficiently large sizes. Important aspects of implementation in distributed architectures are also discussed using the suitability of block SRHT for these architectures.

Block sub-sampled randomized Hadamard transform

For n being a power of two, an SRHT matrix can be defined as follows:

Ω = n l RHD, (2.1)
where R is a l × n uniform, with or without replacement, random sampling matrix, H is a n × n Walsh-Hadamard matrix rescaled by 1 √ n , and D is a diagonal matrix with i.i.d. Rademacher random variables ±1 on the diagonal. The properties of SRHT were thoroughly described in [START_REF] Tropp | Improved analysis of the subsampled randomized Hadamard transform[END_REF] with a follow up analysis in [START_REF] Boutsidis | Improved matrix algorithms via the subsampled randomized hadamard transform[END_REF]. SRHT matrices are commonly used in randomized algorithms as they can be applied to vectors using only n log 2 n flops, while general unstructured matrices require 2nl flops. At the same time, they satisfy the (ε, δ, d) OSE property, if [START_REF] Balabanov | Randomized linear algebra for model reduction. Part I: Galerkin methods and error estimation[END_REF]

l ≥ 3ε -2 ( √ d + 8 log 6n δ ) 2 3d δ , (2.2)
which is only by a logarithmic factor in δ and n larger than the optimal bound. For a general n, a partial SRHT can be used that is defined as the first n columns of an SRHT matrix. Using a partial SRHT is equivalent to padding the input data with zeros to make its dimension a power of two. Unfortunately, products with SRHT matrices are not well suited to distributed computing limiting the benefits of SRHT on modern architectures (see e.g. [START_REF] Yang | Implementing randomized matrix algorithms in parallel and distributed environments[END_REF]). This happens majorly due to computing products with H in tensor form

H = 1 √ n 1 1 1 -1 ⊗ 1 1 1 -1 ⊗ . . . ⊗ 1 1 1 -1 ,
requiring cumbersome reduction operator such as a sequence of arrays of butterflies, rather than a simple addition, which we have with Gaussian matrices. This article attempts to alleviate this problem by constructing Ω block-wise as follows Ω = [Ω (1) Ω (2) . . .

Ω (p) ], (2.3)
where

Ω (i) = r l D (i) RHD (i)
are l × r SRHT matrices related to a unique sampling matrix R and different (independent from each other) diagonal matrices D (i) with i.i.d. Rademacher entries ±1, multiplied from the left by another diagonal matrices D (i) with Rademacher entries, r = n p , 1 ≤ i ≤ p. As in the standard SRHT, the condition that r is a power of two can be achieved by zero padding of the input data. The advantage of Ω defined by (2.3) is that it can be multiplied by an n × d matrix V distributed between p processors with rowwise partitioning, as ΩV =

1≤i≤p

Ω (i) V (i) , (2.4)
where V (i) are the corresponding local blocks of rows of V. In this way, to obtain ΩV one can compute the local contributions Ω (i) V (i) on each processor and then sum-reduce them to the master processor. This makes block SRHT matrices have the same application cost in terms of communication as Gaussian matrices. Thus, they should yield much better scalability of computations than standard SRHT [START_REF] Yang | Implementing randomized matrix algorithms in parallel and distributed environments[END_REF]. The sum-reduce operation requires exchanging O(log p) messages and O(dl log p) per-processor communication volume that can be by a factor O( r l ) less than the volume of communication used by standard SRHT (if l ≤ r). At the same time, block SRHT require less flops per processor than Gaussian matrices. To be more specific, the application cost of block SRHT using (2.4) is only O(rd log r + dl log p) flops per processor, while Gaussian matrices require O(rdl + dl log p) flops per processor. It is deduced that block SRHT matrices are both well-suited to distributed computing and efficient in terms of flops. They are expected to outperform all existing oblivious embedding when the local dimension r and the sampling dimension l are large enough.

The procedure for application of the block SRHT can be easily extended to the case when V is distributed with a 2D partitioning. Namely, to multiply Ω by n × n matrix V distributed over a grid of p × p processors, we first compute the local contributions X (i,j) = Ω (i) V (i,j) on each processor, and then sum-reduce the contributions from the j-th column of blocks to the processor (1, j), 1 ≤ j ≤ p. Note that in this case the resulting matrix Y T = ΩA is distributed with rowwise partitioning over processors (1, 1), (1, 2), . . . , (1, p). This provides the ability to efficiently compute the sketch ΩY, or to orthogonalize Y with a routine suited for distributed computing, with no need to reorganize Y. This can be particularly handy in the low-rank approximation algorithm from Section 3.

We will assume that R in (2.3) samples rows uniformly at random and with replacement. Interestingly, in this case the block SRHT can be viewed as a generalization of the SRHT with replacement and the Rademacher embedding, as it reduces to these maps when r = n and r = 1, respectively. Sampling with replacement can be important, for instance, when r is smaller than the dimension of the embedded subspace.

Theorem 2.1 is the main result of the article. It implies the compatibility of the block SRHT with all randomized methods that rely on OSEs, including the methods in Section 3. The estimate of l in Theorem 2.1 is similar to (2.2) for the standard SRHT matrix, and in particular depends only logarithmically on n and δ.

Theorem 2.1 (Main Theorem). Let 0 < ε < 1 and 0 < δ < 1. Let Ω ∈ R l×n be defined by (2.3). If,

n ≥ l ≥ 3.7ε -2 ( √ d + 4 log n δ + 6.3) 2 log 5d δ ,
then Ω is an (ε, δ, d) OSE.

For better presentation the proof of Theorem 2.1 is deferred to the end of the article (see Section 5).

Randomized low-rank approximation

This section addresses the problem of efficient computation of a rank-k approximation of a large matrix A ∈ R m×n with rapidly decaying spectrum. We provide and analyze some common randomized algorithms for this task. They are presented with a projection-based approach relying solely on the OSE property of the sketching matrix, and are compatible with the block SRHT thanks to Theorem 2.1. A particular focus is given to the scenario where A is uniformly distributed over a 2D grid of processors. For simplicity assume that m ≤ n.

It is a well-known fact that the best rank-k approximation of A in terms of the spectral, trace and Frobenius error is given by

[[A]] k := U k Σ k V T k
, where Σ k is a diagonal matrix of k dominant singular values of A, and U k and V k contain the associated left and right singular vectors. In other words, we have

[[A]] k = arg min rank(B)=k ∥A -B∥ ξ , where ξ = 2, * or F. Furthermore, ∥A -[[A]] k ∥ ξ = σ k+1 if ξ = 2, m i=k+1 σ i if ξ = * , or ( m i=k+1 σ 2 i ) 1 2 if ξ = F, where σ 1 ≥ σ 2 ≥ . . . ≥ σ m denote the singular values of A.
Obtaining the best rank-k approximation can be computationally expensive and often becomes the bottleneck of the algorithm. In such case, one has to turn to alternative methods for computing a low-rank approximation, such as the randomized methods described below.

Randomized Singular Value Decomposition

A low-rank approximation of A can be interpreted as reduction of the range of A to a low-dimensional subspace Q capturing the most of A's action. The SVD approximation [[A]] k corresponds to taking Q as range(U k ). A more efficient way is to take Q = range(AΩ T ), where Ω ∈ R l×n is an OSE [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF][START_REF] David | Sketching as a tool for numerical linear algebra[END_REF]. In this case, the optimal approximation is

[[A]] (RSVD) := arg min range(B)⊆Q ∥A -B∥ ξ . (3.1)
Then a rank-k approximation of A can be obtained by a truncated SVD of ) , which leads to the approximation

[[A]] (RSVD
[[A]] (RSVD) k := [[[[A]] (RSVD) ]] k . Notice that Q T [[A]] (RSVD) k = [[Q T A]] k
, where Q is an orthonormal basis for Q. This observation constitutes the core for the RSVD algorithm (see Algorithm 1) for the computation of Algorithm 8], with specifying the way of computing the SVD of a tall and skinny matrix Z T by a QR factorization and the SVD of the R factor (see steps 4 and 5). The computation of Y can be effectively done with the procedure from Section 2 using the block structure of Ω. The QR factorizations in steps 2 and 4 should be performed with TSQR [START_REF] Demmel | Communication-optimal parallel and sequential qr and lu factorizations[END_REF] or other methods having low communication cost. The computational cost of Algorithm 1 is dominated by computing Q T A in step 3 and possibly AΩ T in step 1.

[[A]] (RSVD) k . Algorithm 1 is same as [1,

Algorithm 1 RSVD, based on [1, Algorithm 8]

Require: m × n matrix A, l × n matrix Ω with l ≪ n, the target rank k.

1: Compute Y = AΩ T . 2: Orthogonalize Y with QR factorization, and get Q.

3: Compute matrix Z = Q T A. 4: Obtain a QR factorization PR of Z T 5: Use SVD to compute the best rank-k approximation Ũk Σk ṼT k of R T . 6: Output factorization [[A]] (RSVD) k = (Q Ũk ) Σk (P Ṽk ) T .
Let us now characterize the accuracy of RSVD approximation. It can be measured for instance by the quasi-optimality constant

∥A-[[A]] (RSVD) k ∥ ξ ∥A-[[A]] k ∥ ξ -1. It is first shown that the accuracy of [[A]] (RSVD) k is guaranteed if Q = range(AΩ T ) captures well the range of A, i.e., for some d ≥ k and ε * ≤ 1 2 we have ∥A -[[A]] (RSVD) ∥ 2 F ≤ (1 + ε * )∥A -[[A]] d ∥ 2 F . (3.2)
Then, by the triangle inequality,

∥A -[[A]] (RSVD) k ∥ ξ ≤ ∥A -[[A]] k ∥ ξ + ∥A -[[A]] (RSVD) ∥ ξ , (3.3)
where ξ = 2, * or F, we obtain

∥A -[[A]] (RSVD) k ∥ ξ ≤ ∥A -[[A]] k ∥ ξ + 2.5∥A -[[A]] d ∥ F . (3.4)
This result guarantees, under the condition (3.2), the quasi-optimality of

[[A]] (RSVD)
k with respect to the Frobenius norm. Furthermore, if A has a fast enough singular value decay and d is large enough, so that the tail after

d-th singular value of A, i.e. ( m i=d+1 σ 2 i ) 1 2 = ∥A -[[A]] d ∥ F is small compared to the k + 1-th singular value σ k+1 = ∥A -[[A]] k ∥ 2 then [[A]] (RSVD)
k is almost as accurate as [[A]] k with respect to all three norms. Moreover, increasing d can make the quasi-optimality constant arbitrary close to zero.

It remains to obtain the conditions on Ω such that (3.2) holds with high probability. This can be done for instance with the results from [START_REF] David | Sketching as a tool for numerical linear algebra[END_REF]. Take ε =

∥V T d Ω T Ω(A -[[A]] d ) T ∥ 2 F ≤ ε∥A -[[A]] d ∥ 2 F .
In turn the latter condition is satisfied with probability at least 1 -δ if Ω is an ( ε d , δ N , 1) OSE, where N = 2md + m + d, as shown below. The OSE property of Ω and the union bound argument guarantee that for given N fixed vectors z i , we have

(1 -ε d )∥z i ∥ 2 2 ≤ ∥Ωz i ∥ 2 2 ≤ (1 + ε d )∥z i ∥ 2 2 , for 1 ≤ i ≤ N (3.5)
with probability at least 1 -δ. Take set {z i } composed of the columns of V d denoted by x i , the columns of

(A -[[A]] d )
T denoted by y i , and all the pairs x i + y j and x i -y j . Then the relation (3.5), the parallelogram identity and the fact that x T i y j = 0, imply that

|x T i Ω T Ωy j | ≤ ε d ∥x i ∥ 2 ∥y j ∥ 2 , for 1 ≤ i ≤ d, 1 ≤ j ≤ m. Consequently, we have ∥V T d Ω T Ω(A -[[A]] d ) T ∥ 2 F = d i=1 m j=1 |x T i Ω T Ωyj| 2 ≤ ε d d i=1 m j=1 ∥xi∥ 2 2 ∥yj∥ 2 2 = ε d ∥V d ∥ 2 F ∥A -[[A]] d ∥ 2 F (3.6)
with probability at least 1 -δ. ). Whereas for Gaussian matrices the required number of rows to satisfy the aforementioned OSEs properties is somewhat lower:

l = O(d log m δ ).
Although it has to be said that SRHT matrices in practice give similar results as Gaussian matrices [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF]. As can be seen from our experiments, this should also be the case for block SRHT. Moreover, we note that the condition l = O(d log n δ ) for Gaussian matrices is still pessimistic. This overestimation is an artifact due to the use of a general analysis based solely on the OSE property. In reality, a Gaussian Ω should satisfy (3.2) with high probability if it has size l = O(d) with a small constant (say 2 or 4) [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF][START_REF] Tropp | Fixed-rank approximation of a positive-semidefinite matrix from streaming data[END_REF].

Nyström approximation

Although being more efficient than the deterministic SVD, the RSVD still can be computationally heavy, especially on distributed architectures, as it requires two passes over A, and a multiplication of A by Q, which is a large dense matrix. Next we discuss improved algorithm that can circumvent these drawbacks. Assume that A is positive semi-definite matrix.

Notice that it can be computationally beneficial to change the norm ∥ • ∥ ξ in (3.1) to its sketched estimate ∥Ω • Ω T ∥ ξ . The accuracy of such an estimation can be guaranteed thanks to the fact that Ω is an OSE. This leads to Nyström approximation [[A]] (Nyst) given below

[[A]] (Nyst) := arg min range(B)⊆Q ∥Ω(A -B)Ω T ∥ ξ , (3.7)
or in a more usual form [START_REF] Alaoui | Fast randomized kernel ridge regression with statistical guarantees[END_REF][START_REF] Tropp | Fixed-rank approximation of a positive-semidefinite matrix from streaming data[END_REF][START_REF] Gittens | Revisiting the nystrom method for improved large-scale machine learning[END_REF][START_REF] Chiu | Sublinear randomized algorithms for skeleton decompositions[END_REF][START_REF] Drineas | On the nyström method for approximating a gram matrix for improved kernel-based learning[END_REF]:

[[A]] (Nyst) = (ΩA) T (ΩAΩ T ) † (ΩA),
where (ΩAΩ T ) † denotes the pseudo-inverse of ΩAΩ T . Then a rank-k approximation of A can be obtained by an SVD of [[A]] (Nyst) , which leads to the approximation

[[A]] (Nyst) k := [[[[A]] (Nyst) ]] k .
This way of obtaining a rank-k approximation from [[A]] (Nyst) is referred to as the modified fixed-rank Nyström via QR [START_REF] Pourkamali | Improved fixed-rank nyström approximation via qr decomposition: Practical and theoretical aspects[END_REF][START_REF] Pourkamali-Anaraki | Randomized clustered nystrom for large-scale kernel machines[END_REF][START_REF] Tropp | Fixed-rank approximation of a positive-semidefinite matrix from streaming data[END_REF] Algorithm 2 describes a way for computing

[[A]] (Nyst)
k suited for distributed computing under 2D partitioning of A. The matrices Y and ΩY can be computed with the procedure from Section 2 using the block structure of Ω. The QR factorization Z = QR in step 4 can be computed with TSQR. Note that in step 6, instead of computing Ûk as (Y Ṽk ) Σ-1 k we could use Ûk = Q Ũk , which would provide more numerical stability but entail a larger computational cost. Algorithm 2 needs only one pass over the matrix A, and does not involve any high-dimensional operations on A except the computation of the sketch Y = AΩ T , which implies its superiority over the standard SVD as well as randomized SVD [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF][START_REF] Tropp | Fixed-rank approximation of a positive-semidefinite matrix from streaming data[END_REF]. In fact, the dominant computational cost of Algorithm 2 is associated with computing Y and ΩY in steps 1 and 2, when r is sufficiently large.

Algorithm 2 Randomized Nyström approximation

Require: n × n matrix A, l × n matrix Ω with l ≪ n, the target rank k.

1: Compute Y = AΩ T .

2: Obtain a Cholesky factor C of ΩY.

3: Compute Z = YC -1 with backward substitution.

4: Obtain the R factor R of Z (with TSQR or similar algorithm). 5: Use SVD to compute the best rank-k approximation Ũk Σk ṼT k of R.

6: Compute Ûk = (Y Ṽk ) Σ-1 k . 7: Output factorization [[A]] (Nyst) k = Ûk Σ2 k ÛT k .
Remark 3.1. The matrix ΩAΩ T can be rank-deficient, for instance, if A or Ω have lower rank than l, which will cause a problem for obtaining a Cholesky factorization in step 2. In this case, a remedy can be to compute an SVD instead of the Cholesky factorization, and take C as a square root of ΩY in SVD form, that then can be used for the pseudo-inversion in step 3. Another possibility is to make A full-rank by using shifting as in [START_REF] Li | Algorithm 971: An implementation of a randomized algorithm for principal component analysis[END_REF].

Let us now characterize the accuracy of Nyström approximation. Notice the following identity [START_REF] Gittens | The spectral norm error of the naive nystrom extension[END_REF]:

A -[[A]] (Nyst) = (A 1 2 -[[A 1 2 ]] (RSVD) ) T (A 1 2 -[[A 1 2 ]] (RSVD) ).
By combining this observation with the derived earlier results on RSVD with A ← A 1 2 , we obtain that

∥A -[[A]] (Nyst) ∥ * = ∥A 1 2 -[[A 1 2 ]] (RSVD) ∥ 2 F ≤ (1 + ε * )∥A 1 2 -[[A 1 2 ]] d ∥ 2 F = (1 + ε * )∥A -[[A]] d ∥ * (3.8)
holds with probability at least 1 -

2δ if Ω is ( 1 3 , δ, d) OSE and ( ε d , δ N , 1) OSE. It is then noticed that (3.

8) also implies the accuracy of the truncated approximation [[A]]

∥A -[[A]] (Nyst) k ∥ ξ ≤ ∥A -[[A]] k ∥ ξ + 2∥A -[[A]] (Nyst) ∥ ξ , (3.9)
where ξ = 2 or * , so that we have by (3.8),

∥A -[[A]] (Nyst) k ∥ ξ ≤ ∥A -[[A]] k ∥ ξ + 3∥A -[[A]] d ∥ * . ( 3

Single-view approximation of non-psd matrix

The Nyström method is applicable only when A is positive semi-definite. Next, we describe a single-view algorithm that works with general matrices. In recent years several such algorithms have been proposed [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF][START_REF] Upadhyay | Fast and space-optimal low-rank factorization in the streaming model with application in differential privacy[END_REF][START_REF] Tropp | Practical sketching algorithms for low-rank matrix approximation[END_REF][START_REF] Tropp | Streaming low-rank matrix approximation with an application to scientific simulation[END_REF][START_REF] Kannan | Principal component analysis and higher correlations for distributed data[END_REF]. The single-view approximation from [START_REF] Tropp | Streaming low-rank matrix approximation with an application to scientific simulation[END_REF] involves four sketching matrices. Two of them, Ω ∈ R l×n , and Γ ∈ R l×m , are used to construct the approximation subspaces Q = range(AΩ T ) and P = range(A T Γ T ) capturing the actions of A and A T . Whereas the other two, Φ ∈ R s×m , and Ψ ∈ R s×n , with l ≤ s, provide the "core sketch" of A. In our projection-based interpretation, the "core sketch" corresponds to estimation of the norm ∥ • ∥ ξ of the residual by ∥Φ • Ψ T ∥ ξ . Thus, we arrive to the following low-rank approximation of A:

[[A]] (sRSVD) = arg min range(B)⊆Q, range(B T )⊆P ∥Φ(A -B)Ψ T ∥ ξ . (3.11) Notice that [[A]] (sRSVD) is given as Q(ΦQ) † (ΦAΨ T )(ΨP) † P T
, where Q denotes an orthogonal basis for Q and P denotes an orthogonal basis for P . Similarly as in the case of RSVD and Nyström approximations, the rank-k approximation of A can be obtained by truncating

[[A]] (sRSVD) with SVD, which provides [[A]] (sRSVD) k = [[[[A]] (sRSVD) ]] k [24]. Notice that [[A]] (sRSVD) k is now given as Q[[C]] k P T , where C = (ΦQ) † (ΦAΨ T )(ΨP) † .

This observation leads to Algorithm 3 for the computation of [[A]] (RSVD)

k . Again, in step 2 one can use TSQR algorithm for the efficient QR factorization on distributed architectures. As in Nyström method, the computational cost of Algorithm 3 is dominated by the applications of sketching matrices.

Algorithm 3 Single-view RSVD [1, Algorithm 17]

Require: m × n matrix A, matrices Ω, Γ with l rows, and Φ, Ψ with s rows, with l ≤ s ≪ m, the target rank k.

1: Compute X = A T Γ T , Y = AΩ T and Z = ΦAΨ T . 2: Orthogonalize X and Y with a QR factorization and obtain Q and P.

To characterize the accuracy of [[A]] (sRSVD)

k we shall assume that Φ, Ψ satisfy both (ε, δ, l) OSE and (ε, δ n , 1) OSE properties with ε ≤ 2 9 . Then these sketching matrices are ε-embeddings for Q and P , i.e. they satisfy (1.1) taking Ω ← Φ, V ← Q or Ω ← Ψ, V ← P , simultaneously, with probability at least 1 -2δ. Moreover, let Π Q and Π P denote the orthogonal projectors onto Q and P . Then Ψ satisfies the ε-embedding property for every subspace spanned by a row of A -Π Q AΠ P simultaneously with probability at least 1 -δ, and Φ satisfies the ε-embedding property for every subspace spanned by a column of (A -Π Q AΠ P )Ψ T simultaneously with probability at least 1 -δ, so that

∥Φ(A -Π Q AΠ P )Ψ T ∥ F ≤ √ 1 + ε∥(A -Π Q AΠ P )Ψ T ∥ F ≤ (1 + ε)∥A -Π Q AΠ P ∥ F
holds with probability at least 1 -2δ. Thus, we have with probability at least 1 -4δ,

∥A -[[A]] (sRSVD) ∥ F ≤ ∥Π Q AΠ P -[[A]] (sRSVD) ∥ F + ∥A -Π Q AΠ P ∥ F ≤ 1 1 -ε ∥Φ(Π Q AΠ P -[[A]] (sRSVD) )Ψ T ∥ F + ∥A -Π Q AΠ P ∥ F ≤ 1 1 -ε ∥Φ(A -[[A]] (sRSVD) )Ψ T ∥ F + 2.6∥A -Π Q AΠ P ∥ F ≤ 1 1 -ε ∥Φ(A -Π Q AΠ P )Ψ T ∥ F + 2.6∥A -Π Q AΠ P ∥ F ≤ 4.2∥A -Π Q AΠ P ∥ F ≤ 4.2(∥A -Π Q A∥ F + ∥A -AΠ P ∥ F ),
which in turn guarantees the accuracy of [[A]] (sRSVD) if Q and P capture well the actions of A and A T . Clearly, ) and ) is the RSVD approximation of A T , associated with the sketching matrix Γ. This suggests that the single-view approximation (3.11) should have a similar quality as the RSVD approximation (3.1). By combining the above consideration with the results on RSVD from Section 3.1 and the triangle inequality we obtain that

Π Q A = [[A]] (RSVD
AΠ P = ([[A T ]] (RSVD) ) T , where [[A T ]] (RSVD
∥A -[[A]] (sRSVD) k ∥ ξ ≤ ∥A -[[A]] (sRSVD) k ∥ ξ + 2∥A -[[A]] (sRSVD) ∥F ≤ ∥A -[[A]] (sRSVD) k ∥ ξ + 16.8 √ 1 + 2.25ε∥A -[[A]] d ∥F
holds with probability at least 1 -6δ for ξ = 2, * or F, given that Φ, Ψ satisfy both (ε, δ, l) OSE and (ε, δ n , 1) OSE properties, and Ω, Γ satisfy both (

Numerical experiments

For numerical experiments, we chose the Nyström approximation as a representative application. The validation of block SRHT is done through comparison with Gaussian embeddings. In the plots, BSRHT refers to block SRHT. The comparison with standard SRHT is impertinent since SRHT matrices are not that well scalable as Gaussian matrices and have no better accuracy [START_REF] Yang | Implementing randomized matrix algorithms in parallel and distributed environments[END_REF].

Nyström approximation

This experiment was executed with Julia programming language version 1.7.2 along with the Distributed.jl and DistributedArrays.jl packages for parallelism. We used 2 nodes Intel Skylake 2.7GHz (AVX512) having 48 available cores and 180 MB of RAM each. In this experiment we used only 32 cores on each node. As input data we took the MNIST or YearPredictionMSD datasets [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF][START_REF] Thierry Bertin-Mahieux | The million song dataset[END_REF]. The radial basis function e -∥xi-xj ∥ 2 /σ 2 was used to build a dense positive definite matrix A of size n × n from n rows of the input data. The parameter σ was chosen as 100 for the MNIST dataset and 10 4 as well as 10 5 for the YearPredictionMSD dataset. The dimension n was taken as 65536.

The matrix A has been uniformly distributed on a square grid of 8 × 8 processors. In all the experiments, the local matrices Ω (i) on each processor were generated with a seeded random number generator with a low communication cost. Figure 1 depicts the convergence of the error of the low-rank approximation obtained with Algorithm 2 taking Ω as a block SRHT. The results for Gaussian Ω are practically identical and therefore are not displayed.

In this numerical experiment, the error is measured with the trace norm. Different sketching sizes l were tested. For each pair of parameters (l, k) 20 different approximations were computed for each type of Ω, in order to have the 95% confidence interval. Nevertheless this interval is not displayed as it is too small to be visible. Figure 2 gives runtime characterization. In particular we depict the runtime spent on computing Y = AΩ T and ΩY in steps 1 and 2 of Algorithm 2. These operations will dominate the overall computational cost, when the block size is large enough. Nevertheless the reader should be aware that TSQR and the SVD of R (step 4 and 5) are also important, especially when the sampling size is close to the block size. The parameter k is not involved in steps 1 and 3 hence not mentioned in Figure 2. According to Figure 2, the runtime of the Gaussian sampling is up to 2.5 times higher and grows faster with l than the runtime taken by block SRHT. Note that for block SRHT the local computation cost is independent of l, hence the slope comes only from the reductions in steps 1 and 2. On the other hand, the Gaussian sampling involves local computations with linear dependency in l, in addition to these reductions.

Cost of application to tall-and-skinny matrix

Next we investigate the performance of block SRHT on larger scale. For this we consider a product of Ω with a tall-and-skinny matrix V, for instance in the context of solving an overdetermined least-squares problem. The same computing environment is used as in the previous experiment, involving now up to 32 nodes and using C99/MPI instead of Julia. The code was compiled using IntelMPI C compiler version 20.0.2 and sequential MKL 20.0.2 with option ILP64. The library FFTW3 used has Intel-specific routines. There is therefore up to 1536 cores available. In this way we generated a random matrix V with d = 200 columns and a variable number n of rows. This matrix was distributed among a variable number p of processors with block rowwise partitioning. Then V was multiplied by either a Gaussian or block SRHT matrix Ω with l = 2000 rows using (2.4). In all experiments, the local Ω (i) matrices on each processor were generated with a seeded random number generator with negligible communication cost. Figure 3 presents a strong scalability test for n = 10 7 . We see that the block SRHT provides an overall speedup by a factor of more than 2.5 over the Gaussian matrices, while demonstrating as good scalability when p ≤ 384. For larger p, however, the reduction operation starts to dominate, which reduces the gain in efficiency. We observe a variability in the MPI_Allreduce operation on larger number of processors for both Gaussian and block SRHT algorithms. However the compute times for both algorithms scale well when increasing the number of processors up to p = 1536. Figure 4 shows a strong scalability test for a higher dimension n = 10 8 . Again we see total" and "BSRHT total" correspond to the overall runtimes, whereas "Gauss. local" and "BSRHT local" stand for the max per-processor runtimes taken by local multiplications. total" and "BSRHT total" correspond to the overall runtimes, whereas "Gauss. local" and "BSRHT local" stand for the max per-processor runtimes taken by local multiplications.

a great scalability of block SRHT for p ≤ 384. For Gaussian matrices, on the other hand, we revealed issues with reaching the memory limit needed to store Ω (i) which made its application on p ≤ 384 processors infeasible. In principle, this problem can be overcome by generating Ω (i) blockwise and applying the blocks to V (i) "on the fly". "Gauss. total" and "BSRHT total" correspond to the overall runtimes, whereas "Gauss. local" and "BSRHT local" stand for the max per-processor runtimes taken by local multiplications.

This however entails a dramatic increase in runtime and therefore is omitted in comparison. On the other hand, for block SRHT we do not have any memory problems 1 . To quantify the advantage of block SRHT in such context, in Figure 5 we provide the memory consumption of the Gaussian and block SRHT matrix. We see that in this sense the reduction in computational cost is indeed drastic. Finally, Figure 6 provides a weak scalability test using n = rp where r = 10 5 . We again see a reduction in runtime of about 2.5 and good scalability for block SRHT up to using p = 1536 processors, similar as in the strong scalability test.

Proof of the main theorem

Before providing the proof for Theorem 2.1, let us first motivate the chosen proof path. Le V be a fixed n × d matrix with orthonormal columns, partitioned using block rowwise partitioning with p blocks V (i) of size r × d.

The statement of the theorem can then be proven by showing that the singular values of ΩV belong to the interval

[ √ 1 -ε, √ 1 + ε] with probability at least 1 -δ. Assume for a moment, that R in (2.
3) is a uniform sampling matrix without replacement. Notice that the random sampling of rows without replacement and then flipping their signs is equivalent to first flipping the signs 1 To reduce the memory consumption, local matrices V (i) are multiplied by Ω (i) in blocks of 20 columns.

where W is a block matrix composed of K × p blocks, with the (j, i)-th block being D (i,j) HD (i) , where D (i,j) are diagonal matrices with entries i.i.d. Rademacher random variables ±1. Unlike WV, the matrix WV (rescaled by 1/ √ K) for sufficiently large K can be proven to be approximately orthonormal with high probability. We are ready to establish the proof of Theorem 2.1.

Notice that the condition in Theorem 2.1 implies that n ≥ l ≥ 3.2ε

(j) = ξ (j) T V/ √ r,
where ξ (j) is a Rademacher vector. Define convex function f

(x) = ∥x T V/ √ r∥ 2 . Observe that f (x) satisfies the Lipschitz bound: ∀x, y, |f (x) -f (y)| ≤ ∥x -y∥ 2 ∥V/ √ r∥ 2 = ∥x -y∥ 2 / √ r.
This allows to apply the Rademacher tail bound P f (ξ (j) ) ≥ Ef (ξ (j) ) + t/ √ r ≤ exp(-t 2 /8), ∀t ≥ 0.

(5.5)

Observe that E(f (ξ (j) )) ≤ (E(f (ξ (j) ) 2 )) Proof. Define τ = ε * /30. Notice that K ≥ 10 4 l ≥ 7.87τ -2 (6.9d + log(r/δ * )).

We have, for any x ∈ R d ,

∥ WVx∥ 2 2 = K j=1 ∥ r i=1
D (i,j) HD (i) V (i) x∥ 2 2 .

(5.6)

Denote by d (k,j) a vector with i-th entry equal to the (k, k)-th entry of matrix D (i,j) , 1 ≤ k ≤ r. Denote by Z (k) the matrix with i-th row equal to the k-th row of matrix HD (i) V (i) , 1 ≤ k ≤ r. Notice the following relations: The statement of the lemma follows by a union bound argument.

Theorem 5.5 (Matrix Chernoff tail bounds from [START_REF] Tropp | Improved analysis of the subsampled randomized Hadamard transform[END_REF]). Consider a finite set X of symmetric positive semi-definite matrices of size d × d. Define the constant L := max Xj ∈X λ max (X j ). Let {X i } l i=1 be a uniformly sampled, without replacement, random subset of X and X := holds with probability at least 1 -4δ * under the condition S. Finally by few algebraic operations, we conclude that, given S, the singular values of ΩV belong to [ √ 1 -ε, √ 1 + ε] with probability at least 4δ/5. The proof of the main theorem is finished by reminding that S occurs with probability at least 1 -δ * , the union bound argument and few additional algebraic operations.

Conclusion

The proposed block SRHT can combine the advantages of structured and unstructured matrices, such as low application complexity and suitability for distributed computing. It should outperform all known embeddings in a distributed architecture with not too large number of processors. At the same time it yields the same approximation guarantees as standard SRHT. We have chosen the low-rank approximation problem as a representative application. We revised popular randomized methods for this problem with implementation aspects on distributed architectures, and then presented their quasi-optimality characterizations from a projection-based point of view, compatible with block SRHT. Numerical validation of the methodology showed that the block SRHT in practice provides solutions of the same quality as Gaussian embeddings. Yet, the block SRHT was up to a factor of 2.5 faster to apply. Moreover, even greater gains in runtime are expected for larger problems and sampling dimensions.

  (Nyst) k due to the following consequence of the triangle inequality (see for instance[START_REF] Tropp | Practical sketching algorithms for low-rank matrix approximation[END_REF] Proposition A.6]):

. 10 )

 10 Similarly to RSVD approximation, this relation guarantees the quasi-optimality of [[A]] (Nyst)k with respect to the trace norm. Furthermore,[[A]] (Nyst)k is almost as accurate as [[A]] k with respect to both the trace norm and the spectral norm, if the tail after d-th singular value of A, i.e. n i=d+1 σ i is small compared to the k + 1-th singular value σ k+1 .

3 : 4 : 5 : 6 :

 3456 Compute ΦQ and ΨP. Compute the core matrix C = (ΦQ) † Z(ΨP) † with least-squares solves. Use SVD to compute the best rank-k approximation Ũk Σk ṼT k of C. Output factorization [[A]] (sRSVD) k = (Q Ũk ) Σk (P Ṽk ).

1 3 ,

 3 δ, d) OSE and ( ε d , δ N , 1) OSE properties, with ε ≤ 2 9 and N = 2nd + n + d. According to Theorem 2.1, these conditions can be satisfied by Ω, Γ that are block SRHT matrices of size l = O(d log 2 n δ ), and Φ, Ψ that are block SRHT matrices of size s = O(d log 3 n δ ).

4 200

 4 Dataset year with σ = 10 Dataset year with σ = 10 5

Figure 1 :

 1 Figure 1: Trace error ∥A -[[A]] (Nyst) k

Figure 2 :

 2 Figure 2: Runtimes of computing Y = AΩ T and ΩY in Algorithm 2 for different sampling sizes.

Figure 3 :

 3 Figure 3: Strong scalability runtimes associated with computing ΩV with n = 10 7 and l = 2000, versus p. "Gauss.total" and "BSRHT total" correspond to the overall runtimes, whereas "Gauss. local" and "BSRHT local" stand for the max per-processor runtimes taken by local multiplications.

Figure 4 :

 4 Figure 4: Strong scalability runtimes associated with computing ΩV with n = 10 8 and l = 2000, versus p. "Gauss.total" and "BSRHT total" correspond to the overall runtimes, whereas "Gauss. local" and "BSRHT local" stand for the max per-processor runtimes taken by local multiplications.

Figure 5 :

 5 Figure 5: Max per-processor memory needed for computing ΩV with n = 10 8 and l = 2000, versus p.

Figure 6 :

 6 Figure 6: Weak scalability runtimes associated with computing ΩV with n = 10 5 × p and l = 2000, versus p."Gauss. total" and "BSRHT total" correspond to the overall runtimes, whereas "Gauss. local" and "BSRHT local" stand for the max per-processor runtimes taken by local multiplications.

1 2 = 1 √K 1 √K

 211 ∥V/ √ r∥ F ≤ d/r. The statement of the lemma follows by combining this relation with (5.5) with t = 8 log 1 δ and using the union bound argument. In Proposition 5.3 is proven that WV with high probability has singular values close to 1. Proposition 5.3. Given S. The singular values of WV with probability at least 1 -δ * lie inside the interval [ 1 -ε * /30, 1 + ε * /30].

D

  (i,j) HD (i) V (i) x∥ 2 2 = r k=1 ⟨d (k,j) , Z (k) x⟩ 2 , 1 ≤ j ≤ K,(5.7)By applying Theorem 5.5 and some algebraic operations, we obtainP(λmin(X)≤ (1 -ε * )λmin(M T M) l rK ) ≤ d e -ε * (1-ε * ) 1-ε * λ min (M T M)l/M ≤ d e -(ε * 2 /2-ε * 3 /6)(M N ) -1 l ≤ δ, P(λmax(X)≥ (1 + ε * )λmax(M T M) l rK ) ≤ d e ε * (1+ε * ) 1+ε * λmax(M T M)l/M ≤ d e -(ε * 2 /2-ε * 3 /6)(M N ) -1 l ≤ δ.

1 √

 1 l i=1 X i . Then P (λ min (X) ≤ (1 -ε)µ min ) ≤ d e -ε (1-ε) 1-ε µmin/L P (λ max (X) ≥ (1 + ε)µ max ) ≤ d e ε (1+ε) 1+ε µmax/Lwhere µ min = lλ min (EX 1 ) and µ max = lλ max (EX 1 ). By plugging Proposition 5.2 and the result of Proposition 5.3 into Proposition 5.4 and taking M = K WV, M = ( √ d + 8 log rK/δ * ) 2 , N = 1.07, along with the union bound argument, we deduce that, √ 1 -ε * 1 -ε * /30 ≤ σ min (ΩV) ≤ σ max (ΩV) ≤ √ 1 + ε * 1 + ε * /30

  4 9 ε * . It follows from [Lemma 45][2] and its proof that (3.2) holds with probability at least 1 -δ if Ω is an ( 1 3 , δ, d) OSE, and

  In Proposition 5.2 is shown that, given S, the matrix WV has rows with equilibrated norms. Notice that, each row of W has entries that are i.i.d Rademacher random variables rescaled by 1/

			√	d + 8 log(rK/δ * )) 2 log(d/δ * ).	(5.4)
	Proposition 5.2. Given S. The rows φ (j) of WV satisfy
	P	max j=1,...,rK	∥φ (j) ∥ 2 ≤ d r + 8 log(rK/δ * )
				√	r.
	Consequently, we have
			φ

* -2 ( r ≥ 1 -δ * .

Proof.

Acknowledgments

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 810367).

and then sampling. By using this consideration, the expression (2.4) can be developed further as

where D (i) are r × r diagonal matrices with Rademacher random variables ±1 on the diagonal, and W = D (1) HD (1) , D (2) HD (2) , . . . , D (l) HD (l) . Looking at (2.4), one can detect many similarities of Ω with standard SRHT matrix. Consequently, in order to argue that ΩV is approximately orthonormal, the first thing to try should be to follow the steps from [START_REF] Tropp | Improved analysis of the subsampled randomized Hadamard transform[END_REF] in the analysis of the original SRHT. In this case the proof recipe would be as follows. First, it could be shown that the matrix W with high probability homogenizes the rows of V. This result then would allow the Matrix Chernoff concentration inequality from [START_REF] Tropp | Improved analysis of the subsampled randomized Hadamard transform[END_REF] to be applied to show that WV and r l RWV have approximately equal minimal and maximal singular values. With these results, it would remain to show that with high probability WV is approximately orthonormal. This, however, can be cumbersome or even impossible in some situations. Think, for example, of the situation when r < d. Therefore, we will assume that R is a uniform sampling matrix with replacement and use the following trick. For better presentation define parameters ε * = 15 16 ε and δ * = δ 5 . Recall that the sampling matrix R restricts a vector x = (x 1 , . . . , x r ) to l coordinates, i.e., we have

The (multi-)set of indices {i 1 , . . . , i l } is a uniform random sample of {1, . . . , r} with replacement. Notice that such sampling of indices is equivalent to the sampling uniformly at random with replacement from {1, . . . , 1, 2, . . . 2, . . . , r, . . . , r} containing K = ⌈10 4 n 2 rδ * ⌉ copies of each index. This observation implies that the sampling matrix R satisfies the identity RH = R[H H . . . H] T = R H, where R is uniform sampling, with replacement, matrix of size l × rK, and H is a block matrix with K blocks of rows, each being equal to H. For a vector x = (x 1 , . . . , x rK ), matrix R satisfies Rx = (x i1 , . . . , x i l ), with 1 ≤ i 1 , . . . , i l ≤ rK, (

where the indices {i 1 , . . . , i l } are drawn uniformly at random with replacement from {1, . . . , rK}. Let S denote the event when i 1 , . . . , i l in (5.2) are all disjoint indices.

Lemma 5.1. S occurs with probability at least 1 -δ * .

Proof. There are in total (rK) l l! ways to select l elements from a (rK)-element set and rK l = rK(rK-1)...(rK-l+1) l!

ways to select l disjoint elements. Consequently, we have

The goal will be to bound the singular values of ΩV under the condition S. The overall probability of success, then will follow by the union bound argument. Next is assumed that S is occurring. Notice that, in this case, matrix R is equivalent to the matrix that samples the entries uniformly at random and without replacement. Then, using the same arguments as in (5.1), we have the following expression for the product ΩV: (5.8)

We have 1

, where Θ is a K × l rescaled Rademacher matrix. By [START_REF] Balabanov | Randomized linear algebra for model reduction. Part I: Galerkin methods and error estimation[END_REF]Proposition 3.7], Θ is an (τ, δ * /r, d) OSE, which implies that

holds with probability at least 1 -δ * /r. By the summation over k and the union bound argument we conclude that

holds with probability at least 1 -δ * . By straightforward substitution of the expressions (5.7) and (5.8) into (5.9), and using (5.6), we conclude that with probability at least

which is equivalent to the statement of the proposition.

Proposition 5.4 presents a corollary of the Matrix Chernoff inequality from [START_REF] Tropp | Improved analysis of the subsampled randomized Hadamard transform[END_REF], used to show that M = 1 √ K WV and rK l RM = ΩV have approximately equal maximal and minimal singular values.

Proposition 5.4 (Corollary of Theorem 2.2 in [START_REF] Tropp | Improved analysis of the subsampled randomized Hadamard transform[END_REF]). Let M be some rK × d matrix. Let 0 < ε * < 1 and 0 < δ * < 1.

Let m (j) denote the rows of M and let M := rK max j=1,...,rK ∥m (j) ∥ 2 2 and N ≥ σ min (M) -2 . Draw at random a sampling matrix R in (5.2) with l ≥ 2(ε * 2 -ε * 3 /3) -1 M N log(d/δ * ). Given S, then with probability at least 1

Proof. For any symmetric matrix X, let λ min (X) and λ max (X) denote the minimal and the maximal eigenvalues of X. To prove Proposition 5.4 we will use the matrix Chernoff tail bounds from [START_REF] Tropp | Improved analysis of the subsampled randomized Hadamard transform[END_REF] presented in Theorem 5.5. Define X := {m (j) (m (j) ) T } n j=1 . Consider the matrix

where T is a set, with #T = l, of elements of {1, 2, . . . , rK} drawn uniformly and without replacement. The matrix X can be written as X = l i=1 X i , where {X i } l i=1 is a uniformly drawn, without replacement, random subset of X. We have E(X 1 ) = 1 rK M T M. Furthermore, λmax(m (j) (m (j) ) T ) = ∥m (j) ∥ 2 ≤ M rK , 1 ≤ j ≤ rK.