Quentin Defenouillere

Yanhua Yao

John Mouaikel

Abdelkader Namane

Aurélie Galopier

Laurence Decourty

Antonia Doyen

Christophe Malabat

Cosmin Saveanu

Quentin Defenouillère

Alain Jacquier
email: jacquier@pasteur.fr

Micheline Fromont-Racine

Cdc48-associated complex bound to 60S particles is required for the clearance of aberrant translation products

Introduction

Cytoplasmic quality control mechanisms prevent the accumulation of potentially toxic translation products by recognition and elimination of non-conforming mRNAs and their associated nascent peptides stalled on ribosomes. Three quality control pathways are known, each targeting a specific type of aberrant mR-NAs: NonSense Mediated mRNA Decay (NMD) [START_REF] Kervestin | NMD: a multifaceted response to premature translational termination[END_REF], NonStop mRNA Decay (NSD) [START_REF] Frischmeyer | An mRNA surveillance mechanism that eliminates transcripts lacking termination codons[END_REF][START_REF] Van Hoof | Exosome-mediated recognition and degradation of mRNAs lacking a termination codon[END_REF] and NoGo mRNA Decay (NGD) (4), eliminating transcripts with a premature termination codon, lacking a STOP codon or carrying a element causing translational arrest, respectively.

NSD constitutes an important quality control process especially in case of [PSI] + state, a physiological context that promotes stop codon read-through, thus generating extended proteins [START_REF] Wilson M A | Genetic interactions between [PSI+] and nonstop mRNA decay affect phenotypic variation[END_REF]. Furthermore, NSD substrates can also result from transcription termination at cryptic polyadenylation (poly(A)) sites located in the coding region of the mRNA. A global map of poly(A) sites recently generated by direct RNA sequencing revealed that many cryptic poly(A) sites exist in yeast as well as in the human genome [START_REF] Ozsolak | Comprehensive polyadenylation site maps in yeast and human reveal pervasive alternative polyadenylation[END_REF].

It has been shown that ribosome stalling and NGD can be induced by the synthesis of highly positively charged aminoacid sequences, such as a poly-arginine or a poly-lysine stretches. The positive charges of these basic tracts would electrostatically interact with the negative charges within the exit tunnel [START_REF] Lu | Electrostatics in the ribosomal tunnel modulate chain elongation rates[END_REF], resulting in ribosome stalling and translation arrest. Such a polybasic tract can arise from the translation of a poly(A) tail into polylysine, which occurs during translation of a Non-Stop mRNA [START_REF] Ito-Harashima | Translation of the poly(A) tail plays crucial roles in nonstop mRNA surveillance via translation repression and protein destabilization by proteasome in yeast[END_REF]. Hence, the boundary between these two quality control pathways remains unclear [START_REF] Shoemaker | Translation drives mRNA quality control[END_REF]. NSD and NGD would both be activated by a stalled ribosome with an empty A site recognized by Dom34/Hbs1, homologs of the eRF1/eRF3 heterodimer, which triggers ribosome dissociation and peptidyl-tRNA release [START_REF] Shoemaker | Dom34:Hbs1 promotes subunit dissociation and peptidyl-tRNA drop-off to initiate no-go decay[END_REF]. In Saccharomyces cerevisiae, Ski7, another homolog of the release factor eRF3, can also trigger NSD [START_REF] Frischmeyer | An mRNA surveillance mechanism that eliminates transcripts lacking termination codons[END_REF][START_REF] Van Hoof | Exosome-mediated recognition and degradation of mRNAs lacking a termination codon[END_REF]. In the case of NGD, a yet unknown endonuclease cleaves the aberrant mRNA. The 5' cleavage product is degraded upon binding of the SKI complex and the exosome whereas the 3' fragment is degraded by Xrn1 (4). It was recently proposed that, in the case of NSD, an endonucleolytic cleavage also occurs, upstream of the ribosome stalled at the poly(A) tail [START_REF] Tsuboi | Dom34:hbs1 plays a general role in quality-control systems by dissociation of a stalled ribosome at the 3' end of aberrant mRNA[END_REF]. In addition, in both NGD and NSD, the nascent peptide is targeted to the proteasome and rapidly degraded. In the case of ribosome stalling induced by the synthesis of consecutive basic amino acid residues, the E3 ubiquitin ligase Not4 was shown to be required for the rapid degradation of the nascent peptide [START_REF] Dimitrova | Nascent peptide-dependent translation arrest leads to Not4p-mediated protein degradation by the proteasome[END_REF], whereas the E3 ligase Ltn1 was found to be required for the ubiquitination and proteasomal degradation of the aberrant peptide synthesized from a Non-Stop mRNA [START_REF] Bengtson | Role of a ribosome-associated E3 ubiquitin ligase in protein quality control[END_REF].

Therefore, during NGD and NSD, the degradation processes of the aberrant mRNA and of the associated nascent peptide are coupled. The absence of the SKI complex combined with the deletion of LTN1 is synthetic slow growth [START_REF] Bengtson | Role of a ribosome-associated E3 ubiquitin ligase in protein quality control[END_REF] showing that a combination of mutations in both Non-Stop mRNA decay and aberrant nascent peptide degradation is detrimental to the cell.

In many cases, such as in ERAD (Endoplasmic-Reticulum-Associated protein Degradation), the AAA-ATPase Cdc48 recognizes polyubiquitinated proteins and is required for their efficient targeting to the proteasome [START_REF] Stolz | Cdc48: a power machine in protein degradation[END_REF][START_REF] Meyer | Emerging functions of the VCP/p97 AAA-ATPase in the ubiquitin system[END_REF]. The binding of Cdc48 to specific targets is determined by its physical interaction with a large number of co-factors. Here we identified two novel factors, Tae2 and Rqc1, which, together with Ltn1 and Cdc48, form a complex associated to the 60S and required for the degradation of nascent peptides synthesized from aberrant Non-Stop mRNAs. We also show that ubiquitination of the nascent peptides starts on the stalled 80S ribosome and continues on the 60S after ribosomal dissociation.

Reserved for Publication Footnotes

Results

A new set of proteins functionally linked to the SKI complex

To identify additional factors involved in the recognition and degradation of aberrant mRNAs or associated nascent peptides, we performed genetic screens using the Genetic Interaction Mapping strategy [START_REF] Decourty | Linking functionally related genes by sensitive and quantitative characterization of genetic interaction profiles[END_REF] with mutants of the SKI complex (ski2∆, ski3∆ and ski8∆) and with ltn1∆. GIM genome-wide genetic screens enable the identification of functional interactions, whether synthetic lethal, synthetic slow growth or epistatic, by combining a query mutation with a collection of mutant strains; either complete gene deletions [START_REF] Giaever | Functional profiling of the Saccharomyces cerevisiae genome[END_REF] or DamP mutants for essential genes [START_REF] Schuldiner | Exploration of the function and organization of the yeast early secretory pathway through an epistatic miniarray profile[END_REF]. Among the twenty mutants having the strongest growth defect when combined with gene deletions of the SKI complex, we identified four gene deletions that were common in all the screens: ltn1∆, vms1∆, ypl009c∆ (tae2∆) and ydr333c∆ (rqc1∆) (Dataset S1). These candidates were also selected in GIM screens performed with ski7∆ (Fig. S1). Figure 1A illustrates the overlap of genetic interactions between ski2∆ and ski3∆. The genetic interactions between the SKI complex and the E3-Ubiquitin ligase Ltn1 confirmed the synthetic slow-growth phenotype reported earlier [START_REF] Bengtson | Role of a ribosome-associated E3 ubiquitin ligase in protein quality control[END_REF]. Vms1 binds Cdc48 and is involved in diverse cellular functions, such as mitochondrial protein turnover [START_REF] Heo | A stress-responsive system for mitochondrial protein degradation[END_REF], ERassociated protein degradation (ERAD) [START_REF] Tran | A Cdc48p-associated factor modulates endoplasmic reticulum-associated degradation, cell stress, and ubiquitinated protein homeostasis[END_REF], or Cdc13 degradation [START_REF] Baek | The cdc48 protein and its cofactor vms1 are involved in cdc13 protein degradation[END_REF]. The deletion of the dubious ORF YDR048C was marked vms1∆* because its deletion overlaps with VMS1 and thus constitutes an independent mutation for this gene. YPL009C was named TAE2 for Translation Associated Element 2, following the identification of this gene deletion in a screen as conferring sensitivity to cycloheximide, a translation inhibitor [START_REF] Alamgir | Chemical-genetic profile analysis of five inhibitory compounds in yeast[END_REF]. YDR333C has been recently named RQC1 for Ribosome Quality Control 1 [START_REF] Brandman | A Ribosome-Bound Quality Control Complex Triggers Degradation of Nascent Peptides and Signals Translation Stress[END_REF]. Remarkably, the GIM screen performed with ltn1∆ selected the four mutants of the SKI complex (ski2∆, ski3∆, ski8∆ and ski7∆) (Fig. 1B).

In the course of this study, we further investigated the cellular role of the two factors Tae2 and Rqc1, which could explain the functional links with the SKI complex. First, we confirmed the functional interactions between TAE2 and RQC1 and the SKI complex by GIM screens with tae2∆ and rqc1∆. Similarly to the ltn1∆ screen, the ski2∆, ski3∆, ski7∆ and ski8∆ mutants exhibited the strongest growth defects in association with tae2∆ or rqc1∆, demonstrating the reciprocity and the specificity of the functional interactions between TAE2, RQC1, LTN1 and the SKI complex (Fig. 1C, 1D and Fig. S1).

We next validated the functional links between the SKI genes and TAE2 or RQC1 by testing the growth of deletion mutant strains on rich medium plates, including ltn1∆ mutant strains as control. Surprisingly, the growth of single or double mutant strains was not significantly affected in these conditions (Fig. 1E). Bengtson and Joazeiro also observed that the deletion of LTN1 had no obvious phenotype on rich medium plates, but demonstrated the hypersensitivity of the ltn1∆ strain to hygromycin B, an antibiotic that decreases translation fidelity and causes stop codon readthrough [START_REF] Bengtson | Role of a ribosome-associated E3 ubiquitin ligase in protein quality control[END_REF]. This observation prompted us to test the viability of the single and double mutant strains in presence of hygromycin B. The growth of each single mutant strain ltn1∆, ski2∆ or tae2∆ was weakly affected but a stronger growth phenotype was observed for rqc1∆. Interestingly, the double deletion strains were clearly more sensitive to hygromycin B when compared to the corresponding simple mutant strains. The ski2∆ mutant exhibited a synthetic slow growth defect in combination with tae2∆ and ltn1∆ and a synthetic lethal effect in combination with rqc1∆ (Fig. 1E). These growth defects correlated with those observed under competitive growth in the GIM screens, which are performed in the presence of G418 and nourseothricin (Fig. 1A-D

Tae2, Rqc1, Ltn1 and Cdc48 are bound to the same 60S ribosomal subunits

The overlap of genetic interaction patterns between TAE2, RQC1 and LTN1 suggested that these factors could be involved in similar biological processes. Since the deletion of any of these genes caused sensitivity to translation-affecting antibiotics such as hygromycin B and cycloheximide [START_REF] Hillenmeyer | The chemical genomic portrait of yeast: uncovering a phenotype for all genes[END_REF], and because it has been shown that Ltn1 binds 60S ribosomal particles and is involved in translational quality control [START_REF] Bengtson | Role of a ribosome-associated E3 ubiquitin ligase in protein quality control[END_REF], we suspected that Tae2 and Rqc1 could be involved in a translation-associated pathway. We therefore analysed the localization of Tae2 and Rqc1 in a polysome gradient to see if they also co-sedimented with the 60S ribosomal particles. We used strains expressing TAP-tagged versions of Tae2, Rqc1 and Ltn1 to monitor the sedimentation of each of these factors in a polysome gradient (Fig. 2A). We observed that Tae2, Rqc1 and Ltn1 were localized in the same fractions as Nog1, a pre60S-binding specific factor here used as a marker for 60S containing fractions, suggesting that all these factors are associated with 60S ribosomal subunits [START_REF] Saveanu | Sequential protein association with nascent 60S ribosomal particles[END_REF].

To see if Tae2, Rqc1 and Ltn1 were bound together to the same 60S particles, we performed Tandem Affinity Purification coupled with label-free mass spectrometry (LC-MS/MS) using strains expressing respectively the Tae2-TAP and Rqc1-TAP fusions as baits and a wild type strain (BY4741) as a negative control. Statistical analysis of label-free quantitative MS data [START_REF] Cox | Andromeda: a peptide search engine integrated into the MaxQuant environment[END_REF] revealed significant enrichment of proteins of the large ribosomal subunit (RPL) in both Tae2-TAP and Rqc1-TAP purifications, confirming the binding of these two factors to 60S particles (Fig. 2B and2C). In addition, we observed that Tae2, Rqc1 and Ltn1 were enriched specifically in both purifications, indicating that they are present together on the same 60S subunits. Moreover, these purifications revealed a strong enrichment of an additional component, Cdc48, and several of its known partners, Ufd1, Npl4, Vms1, Ubx2, Ubx4 and Shp1.

Affinity purifications using Cdc48-TAP as bait with the same wild type strain as a negative control showed a strong enrichment of known Cdc48 partners (Ufd1, Npl4, Ubx2 and Ubx4) and of Ltn1, Tae2 and Rqc1, which confirmed the presence of Cdc48 in complexes that also contain Tae2, Rqc1 and Ltn1 (Fig. 2D and Dataset S2). In parallel, we examined the localization of Cdc48 in a polysome gradient. We observed that, even if most of the Cdc48 hexamers were localized in lighter fractions of the gradient, a subset of Cdc48 sedimented in the 60S fractions (Fig. 2A). This subset is likely to correspond to Cdc48 physically interacting with 60S particles associated with Tae2, Rqc1 and Ltn1.

Ltn1 and Rqc1 are essential for the recruitment of Cdc48 to ribosomal subunits

Cdc48 (called p97 in mammals) is a highly conserved member of the AAA-ATPase family involved in a large number of cellular functions and forming homohexamers able to interact with diverse partners, that confer specificity toward various processes [START_REF] Stolz | Cdc48: a power machine in protein degradation[END_REF]. One of the main functions of Cdc48 is to recognize polyubiquitinated molecules, extract them from a complex and make these substrates accessible to the proteasome for degradation [START_REF] Meyer | Direct binding of ubiquitin conjugates by the mammalian p97 adaptor complexes, p47 and Ufd1-Npl4[END_REF]. Therefore, the ubiquitination activity of Ltn1 could be an important step for the recruitment of the Cdc48-Ufd1-Npl4 complex to 60S particles associated with Tae2 and Rqc1.

We thus tested if the absence of Ltn1 affected the binding of Cdc48 to the Tae2 and Rqc1 associated complexes. We observed that the deletion of LTN1 resulted in a decrease in the amount of Cdc48 and its partners Ufd1 and Npl4 from Tae2-TAP (Fig. 3A and 3B and Dataset S3) and Rqc1-TAP (Fig. 4A and Dataset S4) associated complexes. In addition, we examined the effect of the absence of Rqc1 in Tae2-TAP purifications and noticed a similar Cdc48 decrease (Fig. 3C andDataset S3).

These data were strengthened by an independent experiment showing that the absence of Ltn1 or, to a somewhat lesser extent, Rqc1 efficiently reduced the amount of Cdc48 found in the 60S fractions of a polysome gradient (Fig. S2).

The amount of Cdc48 was also decreased to some extent from the Rqc1-TAP complex in the absence of Tae2 (Fig. 4B and Dataset S4), an observation also corroborated by the diminution of the amount of Cdc48 sedimenting with the 60S in a polysome gradient when Tae2 is absent (Fig. S2). This decrease was paralleled by a similar diminution of Ltn1 (Fig. 4B) suggesting that the decrease of Cdc48 might be indirect, resulting from the lower stability or recruitment of Ltn1 on the complex. In addition, we observed that the Ltn1-3HA fusion, which localizes with the 60S in a polysome gradient in a wild type, is shifted for a large part toward the lighter fractions in the absence of Tae2 (Fig. 4C). This confirms the importance of Tae2 for the stable association of Ltn1 with the complex.

In summary, the Cdc48-Ufd1-Npl4 complex requires both Ltn1 and Rqc1 to efficiently bind 60S particles, which suggests that ubiquitination by Ltn1 is essential for Cdc48 recruitment. Whether Rqc1 has a role in the ubiquitination process as well or if it has a more structural role for Cdc48 binding remains to be determined.

Tae2, Rqc1 and Cdc48 are involved, together with Ltn1, in Non-Stop protein decay

The new 60S associated complex described above contains Cdc48, Rqc1, Tae2 and the E3 ubiquitin Ligase Ltn1, which has been previously shown to ubiquitinate aberrant nascent peptides synthesized from Non-Stop mRNAs to ensure their proteasomal degradation [START_REF] Bengtson | Role of a ribosome-associated E3 ubiquitin ligase in protein quality control[END_REF] essential for the recruitment of the cytoplasmic exosome and thus aberrant mRNA clearance [START_REF] Frischmeyer | An mRNA surveillance mechanism that eliminates transcripts lacking termination codons[END_REF]. Since Tae2 and Rqc1 are functionally linked with the SKI complex and physically associated with Cdc48 and Ltn1 on the 60S, we tested if Tae2, Rqc1 and Cdc48 were also involved in the degradation of aberrant Non-Stop mRNAs and the associated nascent peptides.

The HIS3 Non-Stop reporter (pAV188) (2, 3) was used to test if the Non-Stop degradation pathway was still active in absence of Rqc1 or Tae2. In a wild type strain, his3-Non-Stop mRNAs are efficiently degraded; therefore cells from a parental his3∆100 strain do not grow on media without histidine (Fig. 5A). In contrast, mutations of the SKI complex hinder Non-Stop mRNA decay activity which stabilizes aberrant Non-Stop mRNAs, thus enabling the growth of cells transformed with pAV188 even in absence of histidine. Van Hoof and colleagues showed that a deletion of LTN1 also allowed pAV188 transformed cells to grow in absence of histidine, revealing that the stabilization of proteins synthesized from Non-Stop mRNAs can be sufficient to restore a functional level of His3 activity [START_REF] Wilson M A | A genomic screen in yeast reveals novel aspects of nonstop mRNA metabolism[END_REF]. We performed drop tests using deletion strains for either TAE2 (Fig. 5A, upper panel) or RQC1 (Fig. 5A, lower panel) and observed that both deletions allowed growth in absence of histidine with pAV188. In addition, combinations of TAE2 or RQC1 deletions with ski2∆ (and ltn1∆ in the case of Tae2) revealed a synergistic effect on growth rate on this medium. These results showed that the absence of Tae2 or Rqc1 caused an accumulation of functional His3-Non-Stop proteins, and that this accumulation was even stronger in combination with the absence of Ski2 or Ltn1.

To address whether Tae2 and Rqc1 are involved in Non-Stop mRNA clearance or specifically in the degradation of aberrant nascent peptides, we used the pAV184 reporter which carries a gene coding for Protein A from an mRNA that lacks a termination codon. The same gene but with a normal termination codon was used as a control (pAV183) [START_REF] Frischmeyer | An mRNA surveillance mechanism that eliminates transcripts lacking termination codons[END_REF]. We performed Western blots and Northern blots with total extracts from mutant cells transformed with pAV184 (Non-Stop Protein A) or pAV183 (Stop Protein A). As previously shown, deletion of SKI complex components caused an accumulation of Non-Stop mRNAs and therefore restored expression of aberrant proteins (Fig. S3A &S3B). The ltn1∆ mutant exhibited an accumulation of aberrant proteins without affecting the levels of aberrant Non-Stop mR-NAs (Fig. S3A) [START_REF] Bengtson | Role of a ribosome-associated E3 ubiquitin ligase in protein quality control[END_REF]. We noticed that in the absence of Tae2, the aberrant Non-Stop Protein A was stabilized compared to the wild type (Fig. 5B upper panel) whereas Stop Protein A showed a similar abundance in all tested mutants (Fig. 5B lower panel). We showed that the absence of Tae2 did not affect the level of aberrant Non-Stop mRNA nor its stability (Fig. S3A andS3D). We did not see any clear effect of the RQC1 deletion on Non-Stop Protein A, however, as for the HIS3 Non-Stop reporter, combined deletions of TAE2 or RQC1 with the deletion of SKI2 revealed a synergistic effect on aberrant proteins accumulation (Fig. 5B and5C). We tested in parallel the effect of CDC48 depletion, and also observed a strong accumulation of Non-Stop Protein A (Fig. 5D).

The absence of Tae2, Rqc1 or Cdc48 did not affect the steadystate levels of the Non-Stop Protein A mRNA (Fig. S3A, S3B andS3C), suggesting that like Ltn1, these three factors are specifically involved in the degradation of nascent peptides synthesized from Non-Stop mRNAs but have no role in Non-Stop mRNA decay.

Aberrant nascent peptides remain bound to the ribosome in absence of Tae2, Rqc1 or Cdc48

To understand how Tae2, Rqc1 and Cdc48 act on the recognition or the degradation of aberrant nascent peptides, we analyzed the sedimentation profile of the Non-Stop peptide in a polysome gradient in the absence of Tae2, Rqc1 or Cdc48 (in this experiment, we used a TAP-Non-Stop reporter derived from the Non-Stop Protein A reporter pAV184; see SI Methods). It has been previously shown that in mutants of the SKI complex or in the absence of Ltn1, aberrant nascent peptides accumulated in the 80S fractions of a polysome gradient, suggesting that although Ltn1 sediments with the 60S fractions, the recognition process of aberrant proteins could start on the stalled ribosome and trigger the dissociation of the two ribosomal subunits [START_REF] Bengtson | Role of a ribosome-associated E3 ubiquitin ligase in protein quality control[END_REF]. Whereas the control TAP-Stop protein was localized at the top of the gradient in the free fractions, the TAP-Non-Stop protein accumulated in the 80S fractions in the absence of Tae2, Rqc1 or Cdc48 as in the absence of Ski2 or Ltn1, (Fig. 6A). Interestingly, in the absence of Ltn1 and to a lesser extent in the absence of Rqc1 or Cdc48, we could observe a subset of aberrant proteins sedimenting with the 60S fractions in addition to the 80S fractions. This is consistent with the hypothesis that aberrant nascent peptides could remain blocked in the 60S exit tunnel even after dissociation of the stalled 80S. This subset of TAP-Non-Stop proteins localized in the 60S fractions could thus represent intermediates of degradation between ribosome dissociation and ubiquitination of the aberrant nascent peptide showed by Bengtson M.H. and Joazeiro C. ( 13). Yet, whether the ubiquitination of aberrant nascent peptides starts before or after ribosome dissociation remains unclear.

Cdc48 is involved in the proteasomal targeting of polyubiquitinated substrates stalled on 60S subunits

No high molecular species of the TAP-Non-Stop protein, which could correspond to ubiquitinated forms, were readily apparent in these single-mutant conditions. Given its documented function described in other pathways [START_REF] Stolz | Cdc48: a power machine in protein degradation[END_REF][START_REF] Meyer | Emerging functions of the VCP/p97 AAA-ATPase in the ubiquitin system[END_REF], the Cdc48-Ufd1-Npl4 complex could be expected to play a role in the recognition and proteasomal targeting of polyubiquitinated nascent peptides. We therefore examined the polysome sedimentation profile of TAP-Non-Stop proteins in absence of Ski2 and Cdc48, in order to stabilize TAP-Non-Stop mRNAs and the polyubiquitinated forms of the TAP-Non-Stop proteins, respectively. Besides its strong accumulation in both the 80S and the 60S fractions (Fig. 6B), we observed that upon depletion of Cdc48, a subset of TAP-Non-Stop proteins migrated as high molecular weight proteins, forming a ladder of increasing size from the 80S fractions to the 60S fractions. Affinity purifications of the TAP-Non-Stop reporter and its analysis by Western blot using anti-ubiquitin antibodies revealed the presence of a smear of polyubiquitinated species upon Cdc48 depletion (Fig. 6C), suggesting that the high molecular ladder formed by the TAP-Non-Stop protein in the 60S fractions indeed corresponds to polyubiquitinated forms of this peptide.

Taken together, these results show that i) ubiquitination of the aberrant nascent peptide might start at the 80S level, since some TAP-Non-Stop protein shifted species of intermediate size are already visible in the 80S fractions, ii) high molecular polyubiquitinated forms of the aberrant peptide accumulate on 60S ribosomal particles in absence of Cdc48 and iii) Cdc48 is essential for proteasomal targeting of these polyubiquitinated aberrant substrates, since a depletion of Cdc48 strongly stabilized polyubiquitinated forms of the TAP-Non-Stop protein.

Discussion

Performing GIM genetic screens we identified two new factors, Tae2 and Rqc1, which similarly to Ltn1, become essential during translation stress when the SKI complex is absent (Fig. 1). Remarkably, biochemical purifications revealed that these three factors, together with Cdc48, are part of a common complex in association with the 60S ribosomal subunit (Fig. 2). We found that these factors are all required for the efficient degradation of aberrant nascent peptides synthesized from mRNA lacking a stop codon (NSPD) (Fig. 5).

During the completion of this manuscript, another study was published [START_REF] Brandman | A Ribosome-Bound Quality Control Complex Triggers Degradation of Nascent Peptides and Signals Translation Stress[END_REF] showing the existence of the same Rqc1-Tae2-Ltn1-Cdc48 complex in association with 60S particles, coined the RQC complex for Ribosome Quality Control.

Submission PDF

characterized the RQC complex as a trigger for degradation of stalled polypeptides and as a signal for translational stress, which is fully consistent with our findings.

The presence of a stalled ribosome that does not have a stop codon in its A site was proposed to be the signal that triggers the degradation of Non-Stop or NoGo mRNAs. Dom34/Hbs1 would then be recruited instead of the canonical termination factor eRF1/eRF3, peptidyl-tRNA hydrolysis would not be activated and the nascent peptide, in the form of a peptidyl-tRNA, would remain within the exit tunnel and the P site when Rli1 drives the dissociation of the stalled ribosome. Yet, the absence of Ltn1 was reported to result in the stabilization of the aberrant peptide on the 80S [START_REF] Bengtson | Role of a ribosome-associated E3 ubiquitin ligase in protein quality control[END_REF] suggesting that the ubiquitination process starts before ribosome dissociation. It has been also shown that Ltn1 is localized on the 60S, an observation strengthened by our observation that Ltn1, Tae2, Rqc1 and Cdc48, are present on the same 60S particles, supporting the existence of a hypothetic transient 60S-peptidyl-tRNA intermediate. Indeed, we could detect a Non-Stop peptide bound to the 60S particle (Fig. 6A) and show that polyubiquitinated forms of this peptide accumulate on the 60S in the absence of Cdc48 (Fig. 6B and6C). This observation, together with the fact that Cdc48 does not associate with the RQC particle in the absence of Ltn1 (Fig. 3B and4A), strongly suggests that the Cdc48-Ufd1-Npl4 complex recognizes polyubiquitinated forms of the aberrant peptide at a late stage of the process. Considering the general role of Cdc48, it might extract the polyubiquitinated peptides from the dissociated 60S and escorts them to the proteasome for degradation.

In addition, we observed that Cdc48 binding to the complex is also impeded by the absence of Rqc1 (Fig. 3C), whereas Ltn1 binding is not affected. Whether this reflects the requirement of Rqc1 for peptide ubiquitination or a more structural role of this factor for Cdc48 binding remains to be determined. In addition, the absence of Tae2 weakly affected Cdc48 binding to the complex; this could be an indirect effect due to the decreased binding efficiency of Ltn1 in this context (Fig. 4B and4C).

Finally, Cdc48 was implicated in the degradation of nonfunctional 60S subunits (25S NRD) by promoting the dissociation of stalled 80S, thus allowing the recruitment of the proteasome to stalled 60S to perform ribosomal proteins proteolysis [START_REF] Fujii | 40S subunit dissociation and proteasomedependent RNA degradation in nonfunctional 25S rRNA decay[END_REF]. These results are consistent with the hypothesis that after dissociation of stalled ribosomes, the 60S particles are not recycled to perform new rounds or translation but targeted for proteasomal degradation, together with the aberrant nascent peptide.

Material and methods

Yeast strains, plasmids and GIM screens

The yeast strains used in this study are listed on Table . S1. They were constructed by high-fidelity PCR transformation. Serial dilution were spotted on rich medium plates with or without hygromycin B (40 µg/mL) and incubated for 2 days at 30°C. Details of plasmid constructions and GIM screens ( 16) can be found in SI Methods.

Polysome gradients and Western blotting

Polysome gradients analyses were conducted as described in SI Methods. Collected fractions were precipitated and separated on 10% acrylamide gel. The proteins were detected by hybridization with the appropriate antibodies (Table. S3).

Tandem Affinity Purifications mass spectrometry and data analysis

Biochemical purifications were performed using strains expressing TAP fusion proteins. A first step purification was done with magnetic beads coupled with IgG (SI Methods). An aliquote of the TEV eluate was separated on acrylamide gels and silver stained. The rest of the eluates was analyzed by mass spectrometry. Briefly, after digestion, peptides were identified on a LTQ-Orbitrap and Raw MS data were analyzed using the MaxQuant software as described in SI Methods

Fig. 1 .

 1 Fig. 1. Tae2, Rqc1 and Ltn1 are functionally linked to the SKI complex.(A)GIM screens with ski2∆ and ski3∆ query strains. Each dot indicates a deletion mutant strain. In each screen, we compared the relative changes in a given mutant levels in combination with the query or a reference mutation after 22 generations of a mixture of double mutants population. Microarray estimation of changes were normalised and the log2 of the ration query/reference was plotted. (B) GIM screens with ltn1∆ query strain. As Fig.1A(C) GIM screens with tae2∆ query strain. As Fig.1A(D) GIM screens with rqc1∆ query strain. As Fig.1A(E) Serial dilution of wild type (BY4741), or deletion mutant strainswere spotted on rich medium plates with or without hygromycin B (40 µg/mL) and incubated for 2 days at 30°C.

Fig. 2 .

 2 Fig. 2. Tae2, Rqc1, Ltn1 and Cdc48 are associated to 60S ribosomal particles.(A) Total cellular extracts were prepared from each strain expressing Tae2-TAP, or Rqc1-TAP or Ltn1-TAP and separated on a sucrose gradient (10/30%). Fractions of the gradient were analysed by Western blot using antibodies against TAP, Cdc48, Nog1 and Rps8. (B) Affinity purification of Tae2-TAP associated complex was subjected to LC-MS/MS identification and quantification. Tae2 associated complex was compared to BY4741 strain. Volcano plot shows the fold change (log2 LFQ cpx/LFQ ref.) on the x axis and the p-value distribution (-log10 p-value) on the y axis for the proteins identified in the TAP purification. Each circle indicates an identified protein (in red, for the most significant enriched candidates in the TAP purification; in green for RPLs; in blue for the RPSs).(C) TAP purification with Rqc1-TAP compared to BY4741 strain, As Fig. 2B.(D) TAP purification with Cdc48-TAP compared to BY4741 strain, As Fig. 2B.

Fig. 3 .

 3 Fig. 3. Cdc48 association with Tae2 depends on Ltn1 and Rqc1.(A) Affinity purification of Tae2-TAP associated complex in WT or in the absence of Ltn1 or in the absence of Rqc1. The associated proteins were separated on a polyacrylamide gel and revealed by silver staining. (B) Affinity purification of Tae2-TAP associated complex in presence or absence of Ltn1. The graph shows the comparison between the intensities (LFQ) for each protein identified by LC-MS/MS (log2 scale). (Colour legend as in Fig. 2).(C) Affinity purification of Tae2-TAP associated complex in presence or absence of Rqc1 as in Fig. 3B.

Fig. 4 .

 4 Fig. 4. Ltn1 association to 60S ribosomal particles is affected by Tae2.(A) Affinity purification of Rqc1-TAP associated complex in presence or absence of Ltn1. The graph shows the comparison of the intensities of each protein identified by LC-MS/MS (log2 scale) as in Fig. 3B. (Legend for the dots as in Fig.2A).(B) Affinity purification of Rqc1-TAP associated complex in presence or absence of Tae2 as in Fig.3B. (C) Cellular extracts from the tagged strains Ltn1-3HA, and Ltn1-3HA tae2∆ were separated on a sucrose gradient (10/30%) as in Fig.2. The tagged proteins were revealed by Western blot as in Fig.2using HA antibody.

Fig. 5 .

 5 Fig. 5. The absence of Tae2, Rqc1, Ltn1 as well as Cdc48 leads to the accumulation of Non-Stop proteins.(A) The absence of Tae2 (upper panel) or Rqc1 (lower panel) enables the production of Non-Stop His3. Wild type (BY4741) or deletion mutant strainswere transformed with the Non-Stop HIS3-containing plasmid (pAV188). Serial dilutions of the strains were spotted on SC-Ura and on SC-His and the plates were incubated for three days at 30°C.(B)The same strains as in (A) were transformed with the pAV184 plasmid containing Non-Stop protA gene under the control of GAL10 promoter and with the pAV183 containing protA gene with a Stop codon as control. Non-diluted and 10 folddiluted total cellular extracts were separated on 10% SDS-PAGE polyacrylamide gel and the protein A was revealed by Western blot using PAP antibody. A loading control of the total protein amounts was done using anti-G6PDH antibody. (C) As Fig.5Bwith the strains related to rqc1∆, ski2∆ and rqc1∆ ski2∆.(D) As Fig.5Bwith the strain PrTetO2:CDC48. The depletion of Cdc48 was obtained by incubation of the strain for 12h in presence of Doxycyclin (10 µg/ml).

Fig. 6 .

 6 Fig. 6. Stabilized Non-Stop proteins accumulate on 80S and 60S particles. (A) Cellular extracts from wild type or deletion mutant strains bearing the TAP Non-Stop reporter were separated on sucrose gradients. The sedimentation profile of the TAP-stop protein expressed in wild type strain is shown as control. The fractions were analysed by western blot using antibodies against TAP. In addition, an antibody against Nog1 was used for the WT in order to mark the 60S position. (B) Polysome gradient of cellular extracts from ski2∆ single mutant and ski2∆ Cdc48-depleted strains that express TAP-Non-Stop construct, detected as in (A). (C) Immunoprecipitation of TAP-Non-Stop reporter in ski2∆ single mutant and ski2∆ Cdc48-depleted strain. The cellular extracts from the strains in B were purified on IgG beads, eluted with TEV protease and the associated proteins were analysed by Western blot using monoclonal anti-Ubiquitin antibodies (clone P4D1).

and Fig. S1).

	273
	274
	275
	276
	277
	278
	279
	280
	281
	282
	283
	284
	285
	286
	287
	288
	289
	290
	291
	292
	293
	294
	295 296 297 298 299	Submission PDF
	300
	301
	302
	303
	304
	305
	306
	307
	308
	309
	310
	311
	312
	313
	314
	315
	316
	317
	318
	319
	320
	321
	322
	323
	324
	325
	326
	327
	328
	329
	330
	331
	332
	333
	334
	335
	336
	337
	338
	339
	340
	Footline Author	PNAS	Issue Date	Volume	Issue Number	3

  . In parallel, components of the SKI complex are

	409
	410
	411
	412
	413
	414
	415
	416
	417
	418
	419
	420
	421
	422
	423
	424
	425
	426
	427
	428
	429
	430
	431
	432
	433
	434
	435
	436
	437
	438
	439
	440
	441
	442
	443
	444
	445
	446
	447
	448
	449
	450
	451
	452
	453
	454
	455
	456
	457
	458
	459
	460
	461
	462
	463
	464
	465
	466
	467
	468
	469
	470
	471
	472
	473
	474
	475
	476

www.pnas.org ------Footline Author

www.pnas.org ------Footline Author

Acknowledgments.

We thank A. Van Hoof, A. Buchberger and G. Dieci for providing plasmids and antibodies against Cdc48 and Rps8 respectively. We are grateful to the proteomics platform of the Pasteur Institute for the availability of the Orbitrap. Q.D. and Y.Y. were supported by a fellowship from the Ministère de l'enseignement supérieur et de la Recherche. This work was supported by the grants from the Agence Nationale de la Recherche ANR-2011-BSV6-011-02 and ANR-08-JCJC-0019-01.


Michele Benzi
email: benzi@mathcs.emory.edu

Bora Uçar
email: bora.ucar@ens-lyon.fr

Preconditioning techniques based on the Birkhoff-von Neumann

Keywords: Preconditioning, parallel computing, doubly stochastic matrix, bipartite graphs, Birkhoff-von Neumann decomposition MCS: 65F35, 65F08, 65F50, 65Y05, 15B51, 05C70

We introduce a class of preconditioners for general sparse matrices based on the Birkhoff-von Neumann decomposition of doubly stochastic matrices. These preconditioners are aimed primarily at solving challenging linear systems with highly unstructured and indefinite coefficient matrices. We present some theoretical results and numerical experiments on linear systems from a variety of applications.

Introduction

We consider the solution of linear systems 𝐴𝑥 = 𝑏 where 𝐴 = [𝑎 𝑖𝑗 ] ∈ R 𝑛×𝑛 , 𝑏 is a given vector and 𝑥 is the unknown vector. Our aim is to develop and investigate preconditioners for Krylov subspace methods for solving such linear systems, where 𝐴 is highly unstructured and indefinite.

For a given matrix 𝐴, we first preprocess it to get a doubly stochastic matrix (whose row and column sums are one). Then using this doubly stochastic matrix, we select some fraction of some of the nonzeros of 𝐴 to be included in the preconditioner. Our main tools are the well-known Birkhoff-von Neumann (BvN) decomposition (this will be discussed in Section 2 for completeness), and a splitting of the input matrix in the form 𝐴 = 𝑀 -𝑁 based on its BvN decomposition. When such a splitting is defined, 𝑀 -1 or solvers for 𝑀 𝑦 = 𝑧 are required. We discuss sufficient conditions when such a splitting is convergent and discuss specialized solvers for 𝑀 𝑦 = 𝑧 when these conditions are met. We discuss how to build preconditioners meeting the sufficiency conditions. In case the preconditioners become restrictive in practice, their LU decomposition can be used. Our motivation is that the preconditioner 𝑀 -1 can be applied to vectors via a number of highly concurrent steps, where the number of steps is controlled by the user. Therefore, the preconditioners (or the splittings) can be advantageous for use in many-core computing systems. In the context of splittings, the application of 𝑁 to vectors also enjoys the same property. These motivations are shared by recent work on ILU preconditioners, where their fine-grained computation [START_REF] Chow | Fine-grained parallel incomplete LU factorization[END_REF] and approximate application [START_REF] Anzt | Iterative sparse triangular solves for preconditioning[END_REF] are investigated for GPU-like systems.

The paper is organized as follows. We first give necessary background (Section 2) on doubly stochastic matrices and the BvN decomposition. We then develop splittings for doubly stochastic matrices (Section 3), where we analyze convergence properties and discuss algorithms to construct the preconditioners. Later in the same Section, we discuss how the preconditioners can be used for arbitrary matrices by some preprocessing. Here our approach results in a generalization of the Birkhoff-von Neumann decomposition for matrices with positive and negative entries where the sum of the absolute values of the entries in any given row or column is one. This generalization could be of interest in other areas. Then, we give experimental results (Section 4) with nonnegative and also arbitrary matrices, and then conclude the paper.

Background and definitions

Here we define several properties of matrices: irreducible, full indecomposable, and doubly stochastic matrices.

An 𝑛 × 𝑛 matrix 𝐴 is reducible if there exists a permutation matrix 𝑃 such that

𝑃 𝐴𝑃 𝑇 = [ 𝐴 1,1 𝐴 1,2 𝑂 𝐴 2,2
] ,

where 𝐴 1,1 is an 𝑟 × 𝑟 submatrix, 𝐴 2,2 is an (𝑛 -𝑟) × (𝑛 -𝑟) submatrix, and 1 ≤ 𝑟 < 𝑛. If such a permutation matrix does not exist, then 𝐴 is irreducible [START_REF] Varga | Matrix Iterative Analysis[END_REF]Ch. 1]. When 𝐴 is reducible, either 𝐴 1,1 or 𝐴 2,2 can be reducible as well, and we can recursively identify their diagonal blocks, until all diagonal blocks are irreducible. That is, we can obtain , where each 𝐴 𝑖,𝑖 is square and irreducible. This block upper triangular form, with square irreducible diagonal blocks is called Frobenius normal form [19, p. 532].

𝑃 𝐴𝑃 𝑇 = ⎡ ⎢ ⎢ ⎢ ⎣ 𝐴 1,1 𝐴 1,2 • • • 𝐴 1,𝑠 0 𝐴 2,2
An 𝑛 × 𝑛 matrix 𝐴 is fully indecomposable if there exists a permutation matrix 𝑄 such that 𝐴𝑄 has a zero-free diagonal and is irreducible [START_REF] Brualdi | Combinatorial Matrix Theory, volume 39 of Encyclopedia of Mathematics and its Applications[END_REF]Chs. 3 and 4]. If 𝐴 is not fully indecomposable, but nonsingular, it can be permuted into the block upper triangular form

𝑃 𝐴𝑄 𝑇 = [ 𝐴 1,1 𝐴 1,2 𝑂 𝐴 2,2
] ,

where each 𝐴 𝑖,𝑖 is fully indecomposable or can be further permuted into the block upper triangular form. If the coefficient matrix of a linear system is not fully indecomposable, the block upper triangular form should be obtained, and only the small systems with the diagonal blocks should be factored for simplicity and efficiency [START_REF] Duff | Direct Methods for Sparse Matrices[END_REF]Ch. 6]. We therefore assume without loss of generality that matrix 𝐴 is fully indecomposable. An 𝑛 × 𝑛 matrix 𝐴 is doubly stochastic if 𝑎 𝑖𝑗 ≥ 0 for all 𝑖, 𝑗 and 𝐴𝑒 = 𝐴 𝑇 𝑒 = 𝑒, where 𝑒 is the vector of all ones. This means that the row sums and column sums are equal to one. If these sums are less than one, then the matrix 𝐴 is doubly substochastic. A doubly stochastic matrix is fully indecomposable or is block diagonal where each block is fully indecomposable. By Birkhoff's Theorem [START_REF] Birkhoff | Tres observaciones sobre el algebra lineal[END_REF], there exist coefficients 𝛼 1 , 𝛼 2 , . . . , 𝛼 𝑘 ∈ (0, 1) with ∑ 𝑘 𝑖=1 𝛼 𝑖 = 1, and permutation matrices 𝑃 1 , 𝑃 2 , . . . , 𝑃 𝑘 such that

𝐴 = 𝛼 1 𝑃 1 + 𝛼 2 𝑃 2 + • • • + 𝛼 𝑘 𝑃 𝑘 . ( 1
)
Such a representation of 𝐴 as a convex combination of permutation matrices is known as a Birkhoff-von Neumann decomposition (BvN); in general, it is not unique. The Marcus-Ree Theorem states that there are BvN decompositions with 𝑘 ≤ 𝑛 2 -2𝑛 + 2 for dense matrices; Brualdi and Gibson [START_REF] Brualdi | Convex polyhedra of doubly stochastic matrices: I. Applications of the permanent function[END_REF] and Brualdi [START_REF] Brualdi | Notes on the Birkhoff algorithm for doubly stochastic matrices[END_REF] show that for a fully indecomposable sparse matrix with 𝜏 nonzeros, we have BvN decompositions with 𝑘 ≤ 𝜏 -2𝑛 + 2.

An 𝑛×𝑛 nonnegative, fully indecomposable matrix 𝐴 can be uniquely scaled with two positive diagonal matrices 𝑅 and 𝐶 such that 𝑅𝐴𝐶 is doubly stochastic [START_REF] Sinkhorn | Concerning nonnegative matrices and doubly stochastic matrices[END_REF].

Splittings of doubly stochastic matrices

Definition and properties

Let 𝑏 ∈ R 𝑛 be given and consider solving the linear system 𝐴𝑥 = 𝑏 where 𝐴 is doubly stochastic. Hereafter we assume that 𝐴 is invertible. After finding a representation of 𝐴 in the form [START_REF] Amestoy | A parallel matrix scaling algorithm[END_REF], pick an integer 𝑟 between 1 and 𝑘 -1 and split 𝐴 as

𝐴 = 𝑀 -𝑁, (2)
where

𝑀 = 𝛼 1 𝑃 1 + • • • + 𝛼 𝑟 𝑃 𝑟 , 𝑁 = -𝛼 𝑟+1 𝑃 𝑟+1 -• • • -𝛼 𝑘 𝑃 𝑘 . ( 3
)
Note that 𝑀 and -𝑁 are doubly substochastic matrices.

Definition 1.

A splitting of the form [START_REF] Anzt | Iterative sparse triangular solves for preconditioning[END_REF] with 𝑀 and 𝑁 given by ( 3) is said to be a doubly substochastic splitting.

Definition 2.

A doubly substochastic splitting 𝐴 = 𝑀 -𝑁 of a doubly stochastic matrix 𝐴 is said to be standard if 𝑀 is invertible. We will call such a splitting an SDS splitting.

In general, it is not easy to guarantee that a given doubly substochastic splitting is standard, except for some trivial situation such as the case 𝑟 = 1, in which case 𝑀 is always invertible. We also have a characterization for invertible 𝑀 when 𝑟 = 2.

Theorem 1. Let 𝑀 = 𝛼 1 𝑃 1 + 𝛼 2 𝑃 2 . Then, 𝑀 is invertible if (i) 𝛼 1 ̸ = 𝛼 2 , or (ii) 𝛼 1 = 𝛼 2
and all the fully indecomposable blocks of 𝑀 have an odd number of rows (and columns). If any such block is of even order, 𝑀 is singular.

Proof. We investigate the two cases separately. Case (i): Without loss of generality assume that 𝛼 1 > 𝛼 2 . We have

𝑀 = 𝛼 1 𝑃 1 + 𝛼 2 𝑃 2 = 𝑃 1 (𝛼 1 𝐼 + 𝛼 2 𝑃 𝑇 1 𝑃 2 ).
The matrix 𝛼 1 𝐼 + 𝛼 2 𝑃 𝑇 1 𝑃 2 is nonsingular. Indeed, its eigenvalues are of the form 𝛼 1 + 𝛼 2 𝜆 𝑗 , where 𝜆 𝑗 is the generic eigenvalue of the (orthogonal, doubly stochastic) matrix 𝑃 𝑇 1 𝑃 2 , and since |𝜆 𝑗 | = 1 for all 𝑗 and 𝛼 1 > 𝛼 2 , it follows that 𝛼 1 + 𝛼 2 𝜆 𝑗 ̸ = 0 for all 𝑗. Thus, 𝑀 is invertible.

Case (ii): This is a consequence of the Perron-Frobenius Theorem. To see this, observe that we need to show that under the stated conditions the sum

𝑃 1 + 𝑃 2 = 𝑃 1 (𝐼 + 𝑃 𝑇 1 𝑃 2 ) is invertible, i.e., 𝜆 = -1 cannot be an eigenvalue of 𝑃 𝑇 1 𝑃 2 .
Since both 𝑃 1 and 𝑃 2 are permutation matrices, 𝑃 𝑇 1 𝑃 2 is also a permutation matrix and the Frobenius normal form 𝑇 = Π(𝑃 𝑇 1 𝑃 2 )Π 𝑇 of 𝑃 𝑇 1 𝑃 2 , i.e., the block triangular matrix

𝑇 = ⎡ ⎢ ⎢ ⎢ ⎣ 𝑇 1,1 𝑇 1,2 • • • 𝑇 1,𝑠 0 𝑇 2,2 • • • 𝑇 2,𝑠 . . . . . . . . . . . . 0 0 • • • 𝑇 𝑠,𝑠 ⎤ ⎥ ⎥ ⎥ ⎦ , ( 4
)
has 𝑇 𝑖,𝑗 = 0 for 𝑖 ̸ = 𝑗. The eigenvalues of 𝑃 𝑇 1 𝑃 2 are just the eigenvalues of the diagonal blocks 𝑇 𝑖,𝑖 of 𝑇 . Note that there may be only one such block, corresponding to the case where 𝑃 𝑇 1 𝑃 2 is irreducible. Each diagonal block 𝑇 𝑖,𝑖 is also a permutation matrix. Thus, each 𝑇 𝑖,𝑖 is doubly stochastic, orthogonal, and irreducible. Any matrix of this kind corresponds to a cyclic permutation and has its eigenvalues on the unit circle. If a block has size 𝑟 > 1, its eigenvalues are the 𝑝th roots of unity, 𝜀 ℎ , ℎ = 0, 1, . . . , 𝑝 -1, with 𝜀 = 𝑒 2𝜋𝑖 𝑝 , e.g., by the Perron-Frobenius Theorem (see [17, page 53]). But 𝜆 = -1 is a 𝑝th root of unity if and only if 𝑝 is even. Since 𝑀 is a scalar multiple of 𝑃 1 + 𝑃 2 = 𝑃 1 (𝐼 + 𝑃 𝑇 1 𝑃 2 ), this concludes the proof.

Note that the fully indecomposable blocks of 𝑀 mentioned in the theorem are just connected components of its bipartite graph.

It is possible to generalize the condition (i) in the previous theorem as follows.

Theorem 2. A sufficient condition for 𝑀 =

∑ 𝑟 𝑖=1 𝛼 𝑖 𝑃 𝑖 to be invertible is that one of the 𝛼 𝑖 with 1 ≤ 𝑖 ≤ 𝑟 be greater than the sum of the remaining ones.

Proof. Indeed, assuming (without loss of generality) that

𝛼 1 > 𝛼 2 + • • • + 𝛼 𝑟 , we have 𝑀 = 𝛼 1 𝑃 1 + 𝛼 2 𝑃 2 + • • • + 𝛼 𝑟 𝑃 𝑟 = 𝑃 1 (𝛼 1 𝐼 + 𝛼 2 𝑃 𝑇 1 𝑃 2 + • • • + 𝛼 𝑟 𝑃 𝑇 1 𝑃 𝑟 ).
This matrix is invertible if and only if the matrix

𝛼 1 𝐼 + ∑ 𝑟 𝑖=2 𝛼 𝑖 𝑃 𝑇 1 𝑃 𝑖 is invert- ible. Observing that the eigenvalues 𝜆 𝑗 of ∑ 𝑟 𝑖=2 𝛼 𝑖 𝑃 𝑇 1 𝑃 𝑖 satisfy |𝜆 𝑗 | ≤ ‖𝛼 2 𝑃 𝑇 1 𝑃 2 + • • • + 𝛼 𝑟 𝑃 𝑇 1 𝑃 𝑟 ‖ 2 ≤ 𝑟 ∑ 𝑖=2 𝛼 𝑖 < 𝛼 1
(where we have used the triangle inequality and the fact that the 2-norm of an orthogonal matrix is one), it follows that 𝑀 is invertible.

Again, this condition is only a sufficient one. It is rather restrictive in practice.

Convergence conditions

Let 𝐴 = 𝑀 -𝑁 be an SDS splitting of 𝐴 and consider the stationary iterative scheme

𝑥 𝑘+1 = 𝐻𝑥 𝑘 + 𝑐, 𝐻 = 𝑀 -1 𝑁, 𝑐 = 𝑀 -1 𝑏 , ( 5
)
where 𝑘 = 0, 1, . . . and 𝑥 0 is arbitrary. As is well known, the scheme (5) converges to the solution of 𝐴𝑥 = 𝑏 for any 𝑥 0 if and only if 𝜌(𝐻) < 1. Hence, we are interested in conditions that guarantee that the spectral radius of the iteration matrix

𝐻 = 𝑀 -1 𝑁 = -(𝛼 1 𝑃 1 + • • • + 𝛼 𝑟 𝑃 𝑟 ) -1 (𝛼 𝑟+1 𝑃 𝑟+1 + • • • + 𝛼 𝑘 𝑃 𝑘 )
is strictly less than one. In general, this problem appears to be difficult. We have a necessary condition (Theorem 3), and a sufficient condition (Theorem 4) which is simple but restrictive. with 𝛽 := ∑ 𝑟 𝑖=1 𝛼 𝑖 ∈ (0, 1). Next, observe that 𝑀 -1 𝑁 𝑒 = 𝜆𝑒 is equivalent to 𝑁 𝑒 = 𝜆𝑀 𝑒 or, since we can assume that 𝜆 ̸ = 0, to 𝑀 𝑒 = 1 𝜆 𝑁 𝑒. Substituting 𝛽𝑒 for 𝑀 𝑒 and (𝛽 -1)𝑒 for 𝑁 𝑒 we find

𝛽𝑒 = 1 𝜆 (𝛽 -1)𝑒, or 𝜆 = 𝛽 -1 𝛽 .
Hence, 𝛽-1 𝛽 is an eigenvalue of 𝐻 = 𝑀 -1 𝑁 corresponding to the eigenvector 𝑒. Since |𝜆| = 1-𝛽 𝛽 , we conclude that 𝜌(𝐻) ≥ 1 for 𝛽 ∈ (0, 1 2 ]. This concludes the proof.

Theorem 4. Suppose that one of the 𝛼 𝑖 appearing in 𝑀 is greater than the sum of all the other 𝛼 𝑖 . Then 𝜌(𝑀 -1 𝑁 ) < 1 and the stationary iterative method [START_REF] Brualdi | Notes on the Birkhoff algorithm for doubly stochastic matrices[END_REF] converges for all 𝑥 0 to the unique solution of 𝐴𝑥 = 𝑏.

Proof. Assuming (without loss of generality) that

𝛼 1 > 𝛼 2 + • • • + 𝛼 𝑘 , (6)
(which, incidentally, ensures that the matrix 𝑀 is invertible) we show below that

‖𝐻‖ 2 = ‖(𝛼 1 𝑃 1 + • • • + 𝛼 𝑟 𝑃 𝑟 ) -1 (𝛼 𝑟+1 𝑃 𝑟+1 + • • • + 𝛼 𝑘 𝑃 𝑘 )‖ 2 < 1. (7)
This, together with the fact that 𝜌(𝐻) ≤ ‖𝐻‖ 2 , ensures convergence. To prove [START_REF] Brualdi | Combinatorial Matrix Theory, volume 39 of Encyclopedia of Mathematics and its Applications[END_REF] we start by observing that

𝑀 = 𝛼 1 𝑃 1 + • • • + 𝛼 𝑟 𝑃 𝑟 = 𝛼 1 𝑃 1 ( 𝐼 + 𝛼 2 𝛼 1 𝑄 2 + • • • + 𝛼 𝑟 𝛼 1 𝑄 𝑟 ) ,
where 𝑄 𝑖 = 𝑃 𝑇 1 𝑃 𝑖 . Thus, we have

𝑀 -1 = (𝛼 1 𝑃 1 + • • • + 𝛼 𝑟 𝑃 𝑟 ) -1 = 1 𝛼 1 (𝐼 -𝐺) -1 𝑃 𝑇 1 ,
where we have defined

𝐺 = - 1 𝛼 1 𝑟 ∑ 𝑖=2 𝛼 𝑖 𝑄 𝑖 .
Next, we observe that ‖𝐺‖ 2 < 1, since

‖𝐺‖ 2 ≤ 1 𝛼 1 𝑟 ∑ 𝑖=2 𝛼 𝑖 ‖𝑄 𝑖 ‖ 2 = 𝛼 2 + • • • + 𝛼 𝑟 𝛼 1 < 1
as a consequence of (6). Hence, the expansion

(𝐼 -𝐺) -1 = ∞ ∑ ℓ=0 𝐺 ℓ
is convergent, and moreover

‖(𝐼 -𝐺) -1 𝑃 𝑇 1 ‖ 2 = ‖(𝐼 -𝐺) -1 ‖ 2 ≤ 1 1 -‖𝐺‖ 2 ≤ 1 1 - ( 𝛼2 𝛼1 + • • • + 𝛼𝑟 𝛼1 ) .
The last inequality follows from the fact that ‖𝐺‖ 2 ≤ ∑ 𝑟 𝑖=2 𝛼𝑖 𝛼1 , as can be seen using the triangle inequality and the fact that ‖𝑄 𝑖 ‖ 2 = 1 for all 𝑖 = 2, . . . , 𝑟.

Hence, we have

‖𝑀 -1 𝑁 ‖ 2 ≤ 1 𝛼 1 1 1 - ( 𝛼2 𝛼1 + • • • + 𝛼𝑟 𝛼1 ) ‖𝛼 𝑟+1 𝑃 𝑟+1 + • • • + 𝛼 𝑘 𝑃 𝑘 ‖ 2 .
Using once again the triangle inequality (applied to the last term on the right of the foregoing expression) we obtain

‖𝑀 -1 𝑁 ‖ 2 ≤ 𝛼 𝑟+1 + • • • + 𝛼 𝑘 𝛼 1 -(𝛼 2 + • • • + 𝛼 𝑟 )
.

Using condition [START_REF] Brualdi | Convex polyhedra of doubly stochastic matrices: I. Applications of the permanent function[END_REF] we immediately obtain

𝛼 𝑟+1 + • • • + 𝛼 𝑘 𝛼 1 -(𝛼 2 + • • • + 𝛼 𝑟 ) < 1 ,
therefore ‖𝑀 -1 𝑁 ‖ 2 < 1 and the iterative scheme ( 5) is convergent.

Note that if condition ( 6) is satisfied, then the value of 𝑟 in (3) can be chosen arbitrarily; that is, the splitting will converge for any choice of 𝑟 between 1 and 𝑘. In particular, the splitting

𝐴 = 𝑀 -𝑁 with 𝑀 = 𝛼 1 𝑃 1 and 𝑁 = - ∑ 𝑘 𝑖=2 𝛼 𝑖 𝑃 𝑖 is convergent.
It is an open question whether adding more terms to 𝑀 (that is, using 𝑀 = ∑ 𝑟 𝑖=1 𝛼 𝑖 𝑃 𝑖 with 𝑟 > 1) will result in a smaller spectral radius of 𝐻 (and thus in faster asymptotic convergence); notice that adding terms to 𝑀 will make application of the preconditioner more expensive.

Note that condition ( 6) is very restrictive. It implies that 𝛼 1 > 1/2, a very strong restriction. It is clear that given a doubly substochastic matrix, in general it will have no representation (1) with 𝛼 1 > 1/2. On the other hand, it is easy to find examples of splittings of doubly substochastic matrices for which 𝛼 1 = 1/2 and the splitting (3) with 𝑟 = 1 is not convergent. Also, we have found examples with 𝑘 = 3, 𝑟 = 2 and 𝛼 1 + 𝛼 2 > 𝛼 3 where the splitting did not converge. It is an open problem to identify other, less restrictive conditions on the 𝛼 𝑖 (with 1 ≤ 𝑖 ≤ 𝑟) that will ensure convergence, where the pattern of the permutations could also be used.

Solving linear systems with 𝑀

The stationary iterations [START_REF] Brualdi | Notes on the Birkhoff algorithm for doubly stochastic matrices[END_REF] for solving 𝐴𝑥 = 𝑏 or Krylov subspace methods using 𝑀 as a preconditioner require solving linear systems of the form 𝑀 𝑧 = 𝑦.

Assume that 𝑀 = 𝛼 1 𝑃 1 + 𝛼 2 𝑃 2 and 𝛼 1 > 𝛼 2 . The stationary iterations 𝑧 𝑘+1 = 1 𝛼1 𝑃 𝑇 1 (𝑦 -𝛼 2 𝑃 2 𝑧 𝑘
) are convergent for any starting point, with the rate of 𝛼2 𝛼1 . Therefore, if 𝑀 = 𝛼 1 𝑃 1 + 𝛼 2 𝑃 2 and 𝑀 is as described in Theorem 1(i), then we use the above iterations to apply 𝑀 -1 (that is, solve linear systems with 𝑀 ). If 𝑀 is as described in Theorem 4, then we can still solve 𝑀 𝑧 = 𝑦 by stationary iterations, where we use the splitting 𝑀 = 𝛼 1 𝑃 1 -∑ 𝑘 𝑟=2 𝛼 𝑟 𝑃 𝑟 . We note that application of ∑ 𝑘 𝑟=2 𝛼 𝑟 𝑃 𝑟 to a vector 𝑦, that is 𝑧 = ( ∑ 𝑘 𝑟=2 𝛼 𝑟 𝑃 𝑟 )𝑦 can be effectively computed in 𝑘 -1 steps, where at each step, we perform 𝑧 ← 𝑧+(𝛼 𝑟 𝑃 𝑟 )𝑦. This last operation takes 𝑛 input entries, scales them and adds them to 𝑛 different positions in 𝑧. As there are no read or write dependencies between these 𝑛 operations, each step is trivially parallel; especially in shared memory environments, the only parallelization overhead is the creation of threads. Either input can be read in order, or the output can be written in order (by using the inverse of the permutation 𝑃 𝑟 ).

As said before, splitting is guaranteed to work when 𝛼 1 > 1/2. Our experience showed that it does not work in general when 𝛼 1 < 1/2. Therefore, we suggest using 𝑀 as a preconditioner in Krylov subspace methods. There are two ways to do that. The first is to solve linear systems with 𝑀 with a direct method by first factorizing 𝑀 . Considering that the number of nonzeros in 𝑀 would be much less than that of 𝐴, factorization of 𝑀 could be much more efficient than factorization of 𝐴. The second alternative, which we elaborate in Section 3.4, is to build 𝑀 in such a way that one of the coefficients is larger than the sum of the rest.

Algorithms for constructing the preconditioners

It is desirable to have a small number 𝑘 in the Birkhoff-von Neumann decomposition while designing the preconditioner. This is because of the fact that if we use splitting, then 𝑘 determines the number of steps in which we compute the matrix vector products. If we do not use all 𝑘 permutation matrices, having a few with large coefficients should help to design the preconditioner. The problem of obtaining a Birkhoff-von Neumann decomposition with the minimum number 𝑘 is a difficult one, as pointed out by Brualdi [START_REF] Brualdi | Notes on the Birkhoff algorithm for doubly stochastic matrices[END_REF] and has been shown to be NP-complete [START_REF] Dufossé | Notes on Birkhoff-von Neumann decomposition of doubly stochastic matrices[END_REF]. This last reference also discusses a heuristic which delivers a small number of permutation matrices. We summarize the heuristic in Fig. 1.

As seen in Fig. 1, the heuristic proceeds step by step, where at each step 𝑗, a bottleneck matching is found to determine 𝛼 𝑗 and 𝑃 𝑗 . A bottleneck matching can be found using MC64 [START_REF] Duff | The design and use of algorithms for permuting large entries to the diagonal of sparse matrices[END_REF][START_REF] Duff | On algorithms for permuting large entries to the diagonal of a sparse matrix[END_REF]. In this heuristic, at Line 4, 𝑀 is a bottleneck perfect matching; that is, the minimum weight of an edge in 𝑀 is the maximum among all minimum elements of perfect matchings. Also, at Line 5, 𝛼 𝑘 is equal to the bottleneck value of the perfect matching 𝑀 . A nice property of this heuristic is that it delivers 𝛼 𝑖 s in a non-increasing order; that is 𝛼 𝑗 ≥ 𝛼 𝑗+1 for Input: 𝐴, a doubly stochastic matrix; Output: a BvN decomposition (1) of 𝐴 with 𝑘 permutation matrices.

(1) 𝑘 ← 0 (2) while (nnz(𝐴) > 0)

(3) 𝑘 ← 𝑘 + 1 (4)
𝑃 𝑘 ← the pattern of a bottleneck perfect matching 𝑀 in 𝐴 [START_REF] Brualdi | Notes on the Birkhoff algorithm for doubly stochastic matrices[END_REF] 𝛼 𝑘 ← min 𝑀 𝑖,𝑃𝑘 (𝑖)

(6) 𝐴 ← 𝐴 -𝛼 𝑘 𝑃 𝑘 (7) endwhile
all 1 ≤ 𝑗 < 𝑘. The worst case running time of a step of this heuristic can be bounded by the worst case running time of a bottleneck matching algorithm. For matrices where nnz(𝐴) = 𝑂(𝑛), the best known algorithm is of time complexity 𝑂(𝑛 √ 𝑛 log 𝑛), see [START_REF] Gabow | Algorithms for two bottleneck optimization problems[END_REF]. We direct the reader to [8, p. 185] for other cases. This heuristic could be used to build an 𝑀 such that it satisfies the sufficiency condition presented in Theorem 4. That is, we can have 𝑀 with 𝛼1 ∑ 𝑘 𝑖=1 𝛼𝑖 > 1/2, and hence 𝑀 -1 can be applied with splitting iterations. For this, we start by initializing 𝑀 to 𝛼 1 𝑃 1 . Then, when the 𝛼 𝑗 where 𝑗 ≥ 2 is obtained at Line 5, we add 𝛼 𝑗 𝑃 𝑗 to 𝑀 if 𝛼 1 is still larger than the sum of the other 𝛼 𝑗 included in 𝑀 . In practice, we iterate the while loop until 𝑘 is around 10 and collect 𝛼 𝑗 's as described above as long as

𝛼1 ∑ 𝑃 𝑗 ∈𝑀 𝛼𝑗 > 1
1.9 ≈ 0.53.

Arbitrary coefficient matrices

Here we discuss how to apply the proposed preconditioning technique to all fully indecomposable sparse matrices.

Let 𝐴 ≥ 0 be an 𝑛 × 𝑛 fully indecomposable matrix. The first step is to scale the input matrix with two positive diagonal matrices 𝑅 and 𝐶 so that 𝑅𝐴𝐶 is doubly stochastic, or nearly so. For this step, there are different approaches [START_REF] Knight | A fast algorithm for matrix balancing[END_REF][START_REF] Knight | A symmetry preserving algorithm for matrix scaling[END_REF][START_REF] Sinkhorn | Concerning nonnegative matrices and doubly stochastic matrices[END_REF]. Next the linear system 𝐴𝑥 = 𝑏 can be solved by solving 𝑅𝐴𝐶𝑥 ′ = 𝑅𝑏 and recovering 𝑥 = 𝐶𝑥 ′ , where we have a doubly stochastic coefficient matrix.

Suppose 𝐴 is an 𝑛 × 𝑛 fully indecomposable matrix with positive and negative entries. Then, let 𝐵 = abs(𝐴) and consider 𝑅 and 𝐶 making 𝑅𝐵𝐶 doubly stochastic, which has a Birkhoff-von Neumann decomposition 𝑅𝐵𝐶 = ∑ 𝑘 𝑖=1 𝛼 𝑖 𝑃 𝑖 . Then, 𝑅𝐴𝐶 can be expressed as

𝑅𝐴𝐶 = 𝑘 ∑ 𝑖=1 𝛼 𝑖 𝑄 𝑖 .
where

𝑄 𝑖 = [𝑞 (𝑖) 𝑗𝑘 ] 𝑛×𝑛 is obtained from 𝑃 𝑖 = [𝑝 (𝑖)
𝑗𝑘 ] 𝑛×𝑛 as follows:

𝑞 (𝑖) 𝑗𝑘 = sign(𝑎 𝑗𝑘 )𝑝 (𝑖) 𝑗𝑘 .
That is, we can use a convex combination of a set of signed permutation matrices to express any fully indecomposable matrix 𝐴 where abs(𝐴) is doubly stochastic. We can then use the same construct to use 𝑟 permutation matrices to define 𝑀 (for splitting or for defining the preconditioner). We note that the Theorems 2 and 4 remain valid without changes, since the only property we use is the orthogonality of the 𝑃 𝑖 (not the nonnegativity). Theorem 1, on the other hand, needs some changes. All we can prove in this more general setting is that if

𝛼 1 ̸ = ±𝛼 2 , then 𝑀 = 𝛼 1 𝑃 1 + 𝛼 2 𝑃 2 is nonsingular.
We could not find a simple condition for the case 𝛼 1 = ±𝛼 2 since we cannot use the Perron-Frobenius Theorem to conclude anything about 1 (or -1) being an eigenvalue of 𝑃 𝑇 1 𝑃 2 .

Experiments

We are interested in testing our preconditioner on challenging linear systems that pose difficulties for standard preconditioned Krylov subspace methods. We conducted experiments with the preconditioned GMRES of Matlab, without restart. We also use an implementation of the flexible GMRES (FGMRES) [START_REF] Saad | A flexible inner-outer preconditioned GMRES algorithm[END_REF] when needed. We used a large set of matrices which come from three sources. The first set contains 22 matrices which were used by Benzi et al. [START_REF] Benzi | Preconditioning highly indefinite and nonsymmetric matrices[END_REF]. Another six matrices were used by Manguoglu et al. [START_REF] Manguoglu | Weighted matrix ordering and parallel banded preconditioners for iterative linear system solvers[END_REF]; they have experiments with larger matrices but we chose only those with 𝑛 ≤ 20000. These matrices are shown in Table 1. All these matrices, except slide, two-dom, watson4a, and watson5a, are available from the UFL Sparse Matrix Collection [START_REF] Davis | The University of Florida sparse matrix collection[END_REF]. To this set, we add all real, square matrices from the UFL Sparse Matrix Collection which contain "chemical" as the keyword. These matrices pose challenges to Krylov subspace methods. There were a total of 70 such matrices; taking the union with the previously described 28 matrices yield 90 matrices in total. Table 1 shows the size and the number of nonzeros of the largest fully indecomposable block (these are the largest square blocks in the Dulmage-Mendelsohn decomposition of the original matrices) of the first two sets of matrices. The experiments are conducted with those largest blocks; from now on a matrix refers to its largest block.

In Sections 4.1 and 4.2, we present two sets of experiments: with nonnegative matrices and with general matrices. Each time, we compare three preconditioners. The first is ILU(0), which forms a basis of comparison. We created ILU(0) preconditioners by first permuting the coefficient matrices to have a maximum product matching in the diagonal (using MC64), as suggested by Benzi et al. [START_REF] Benzi | Preconditioning highly indefinite and nonsymmetric matrices[END_REF], and then calling ilu of Matlab with the option nofill. The second set of preconditioners is obtained by a (generalized) BvN decomposition and taking a few different values of "𝑟" in defining the preconditioner. We use BvN 𝑟 to denote these preconditioners, where the first 𝑟 permutation matrices and the corresponding coefficients are taken from the output of Algorithm 1. Since we cannot say anything special about these preconditioners, we use their LU decomposition to apply the preconditioner. For this purpose, we used MATLAB's LU with four output arguments (which permutes the matrices for numerical stability and fill-in). The third set of preconditioners is obtained by the method proposed in Section 3.4 so that the preconditioners satisfy the sufficiency condition of Theorem 4, and hence we can use the splitting based specialized solver. We use BvN * to denote these preconditioners. In the experiments with the ILU(0) and BvN 𝑟 preconditioners, we used MATLAB's GMRES, whereas with BvN * we used FGMRES, since we apply BvN * with the specialized splitting based solver. In all cases, we asked a relative residual of 10 -6 from (F)GMRES and run them without restart with at most 3000 iterations. We checked the output of (F)GMRES, and marked those with flag̸ = 0 as unsuccessful. For the successful runs, we checked the relative residual ‖𝐴𝑥-𝑏‖ ‖𝑏‖ for a linear system 𝐴𝑥 = 𝑏, and deemed the runs whose relative residual is larger than 10 -4 as unsuccessful.

Nonnegative matrices

The first set of experiments is conducted on nonnegative matrices. Let 𝐴 be a matrix from the data set (Table 1 and another 62 matrices) and 𝐵 = abs(𝐴), that is, 𝑏 𝑖𝑗 = |𝑎 𝑖𝑗 |. We scaled 𝐵 to a doubly stochastic form with the method of Knight and Ruiz [START_REF] Knight | A fast algorithm for matrix balancing[END_REF]. We asked a tolerance of 10 -8 (so that row and column sums can deviate from 1 by 10 -8 ). We then obtained the Birkhoff-von Neumann decomposition by using Algorithm 1. When the bottleneck value found at a step was smaller than 10 -10 , we stopped the decomposition process-hence we obtain an "approximate" Birkhoff-von Neumann decomposition of an "approximately" doubly stochastic matrix. This way we obtained Then, let 𝑥 ⋆ be a random vector whose entries are from the uniform distribution on the open interval (0, 1), generated using rand of Matlab. We then defined 𝑏 = 𝐵𝑥 ⋆ to be the right hand side of 𝐵𝑥 = 𝑏. We report our experiments with GMRES and FGMRES using different permutation matrices on nonnegative matrices in Table 2. The first part of the table gives the number of GMRES iterations with the ILU(0) preconditioner, and with the proposed Birkhoff-von Neumann based preconditioner, using 𝑟 = 1, 2, 4, 8, 16, 32, and 64 permutation matrices and their coefficients. When there were less than 𝑟 permutation matrices in the BvN decomposition, we used all available permutation matrices. Note that this does not mean that we have an exact inverse, as the BvN decomposition is only approximative. We also report results with FGMRES for the BvN * preconditioners. For this case, we give the number 𝑟 of permutation matrices used and the number of FGMRES iterations (under the column "it"). In the second part (the last row of the table), we give the number of successful runs with different preconditioners; here we also give the average number of permutation matrices in BvN * under the column "𝑟". Some observations are in order. As seen in Table 2, ILU(0) results in 11 successful runs on the matrices listed in the first part of the table, and 32 successful runs in all 90 matrices. BvN 𝑟 preconditioners with differing number of permutation matrices result in at least 53 successful runs in total, where BvN * obtained the highest number of successful runs, but it has in only a few cases the least number of iterations. In 20 cases in the first part of the table, we see that adding more permutation matrices usually helps in reducing the number of iterations. This is not always the case though. We think that this is due to the preconditioner becoming badly conditioned with the increasing number of permutation matrices; if we use the full BvN decomposition, we will have the same conditioning as in 𝐴.

𝐵 ≈ 𝛼 1 𝑃 1 + 𝛼 2 𝑃 2 + • • • + 𝛼 𝑘 𝑃 𝑘 .
In order to clarify more, we present the performance profiles [START_REF] Dolan | Benchmarking optimization software with performance profiles[END_REF] for ILU(0), BvN 𝑟 with 𝑟 = 2, 4, 8, and BvN * in Fig. 2a. A performance profile for a preconditioner shows the fraction of the test cases in which the number of (F)GMRES iterations with the preconditioner is within 𝜏 times the smallest number of (F)GMRES iterations observed (with the mentioned set of preconditioners). Therefore, the higher the profile of a preconditioner, the better is its performance. As seen in this figure, BvN * and BvN 8 manifest themselves as the best preconditioners; they are better than others after around 𝜏 = 1.1, while BvN 8 being always better than others. BvN * is a little behind ILU(0) at the beginning but then catches up with BvN 8 at around 1.2 and finishes as the best alternative.

Although we are mostly concerned with robustness of the proposed preconditioners and their potential for parallel computing, we give a few running time Table 2. The number of GMRES iterations with ILU(0) and BvN𝑟 with different number of permutation matrices, and the number of permutation matrices and FGMRES iterations with BvN*. (F)GMRES are run with tolerance 10 -6 , without restart and with at most 3000 iterations. The symbol "-" flags the cases where (F)GMRES were unsuccessful. All matrices are nonnegative.

BvN𝑟 (iterations)

BvN results on a sequential Matlab environment (on a core of an Intel Xeon E5-2695 with 2.30 GHz clock speed). We present the complexity of the preconditioners and the running time of (F)GMRES (total time spent in the iterations) for the cases where ILU(0) resulted in convergence in Table 3. This set of matrices is chosen to give the running time of (F)GMRES with all preconditioners under study. Otherwise, it is not very meaningful to use sophisticated preconditioners when ILU(0) is effective. Furthermore, the applications of ILU(0) and BvN 𝑟 require triangular solves, which are efficiently implemented in MATLAB. On the other hand, the application of BvN * requires permutations and scaling, which should be very efficient in parallel. The running time given in the right side of Table 3 should not be taken at face value. The complexities of BvN 𝑟 are given as the ratio (nnz(𝐿+𝑈 )-𝑛)/ nnz(𝐴), where 𝐿 and 𝑈 are the factors of the preconditioner. In this setting, ILU(0) has always a complexity of 1.0. The complexity of BvN * is given as nnz(𝑀 )/ nnz(𝐴). As seen in this table, the LU-factors of BvN 𝑟 with 𝑟 ≤ 8 have reasonable number of nonzeros; the complexity of the preconditioners is less than one in all cases except two-dom, lns_3937, west0655, and appu. In some preliminary experiments, we have seen large numbers when 𝑟 ≥ 16. This was especially important for appu where the full LU factorization of 𝐴 contained 85.67 • nnz(𝐴) nonzeros and 𝑀 suffered large fill-in. We run the solver with a large error tolerance of 1.0e-1, as this was enough to get the outer FGMRES to converge. The outer FGMRES iterations would change if a lower error tolerance is used in the inner solver, but we do not dwell into this issue here. With the specified parameters, an application of BvN * required, for the data shown in Table 3, between 128 and 162 iterations.

Table 3. The complexity of the preconditioners and the running time of the solver. The complexities of BvN𝑟 are given as the ratio (nnz(𝐿 + 𝑈 ) -𝑛)/ nnz, where 𝐿 and 𝑈 are the factors of the preconditioner (ILU(0) has a complexity of 1.0), and that of BvN* is given as nnz(𝑀 )/ nnz(𝐴). A sign of "-" in the running time column flags the cases where (F)GMRES were unsuccessful. All matrices are nonnegative. We now comment on the running time of the preconditioner set up phase. The first step is to apply the scaling method of Knight and Ruiz [START_REF] Knight | A fast algorithm for matrix balancing[END_REF]. The dominant cost in this step is sparse matrix-vector and sparse matrix transposevector multiply operations. The second step is to obtain a BvN decomposition, or a partial one, using the algorithms for the bottleneck matching problem from MC64 [START_REF] Duff | The design and use of algorithms for permuting large entries to the diagonal of sparse matrices[END_REF][START_REF] Duff | On algorithms for permuting large entries to the diagonal of a sparse matrix[END_REF]; these algorithms are highly efficient. The number of sparse matrix-vector and sparse matrix transpose-vector multiply operations (MVP) in the scaling algorithm, the running time of the scaling algorithm with the previous setting, and the running time of the partial BvN decomposition with 16 permutation matrices are shown for the same set of matrices in Table 4. Since an application of BvN * required between 128 and 162 iterations, and those iterations are no more costly than sparse matrix vector operations, we could relate the running time of FGMRES with BvN * in Table 3 to the scaling algorithm. For example, for lns_3937 a matrix vector multiplication should not cost more than 0.27/3732 seconds, while the application of BvN * required a total of 1138 • 152 steps of inner solver. In principle the total cost of the inner solver should be around 12.5, but in Table 3, we see 519.83 seconds of total FGMRES time. As in Table 4, for this set of matrices the scaling algorithm and the partial BvN decomposition algorithms are fast. There are efficient parallel scaling algo- rithms [START_REF] Amestoy | A parallel matrix scaling algorithm[END_REF], and the one that we used [START_REF] Knight | A fast algorithm for matrix balancing[END_REF] can also be efficiently parallelized by tapping into the parallel sparse matrix vector multiplication algorithms.

General matrices

In this set of experiments, we used the matrices listed before, while retaining the sign of the nonzero entries. We used the same scaling algorithm as in the previous section and the proposed generalized BvN decomposition (Section 3.5) to construct the preconditioners. As in the previous subsection, we report our experiments with GMRES and FGMRES using different permutation matrices in Table 5. In particular, we report the number of GMRES iterations with the ILU(0) preconditioner, and with the proposed Birkhoff-von Neumann based preconditioner, using 𝑟 = 1, 2, 4, 8, 16, 32, and 64 permutation matrices and their coefficients. We also report results with FGMRES for the BvN * preconditioners. Again, the number 𝑟 of permutation matrices used and the number of FGMRES iterations (under the column "it.") are given for BvN * . The last row of the table gives the number of successful runs with different preconditioners in all 90 matrices; here the average number of permutation matrices in BvN * is also given under the column "𝑟". We again use the performance profiles shown in Fig. 2b to see the effects of the preconditioners. As seen in this table and the associated performance profile, the preconditioners behave much like they do in the nonnegative case. In particular, BvN * and BvN 8 deliver the best performance in terms of number of iterations, followed by ILU(0). We also give the complexity of the preconditioners and the running time of (F)GMRES (total time spent in the iterations) for the cases where ILU(0) resulted in convergence in Table 6. The complexities of the preconditioners remain virtually the same, as expected. As before, the running time is just to give a rough idea in a sequential MATLAB environment. We note that BvN * required, on average, between 125 and 151 inner iterations. Finally note that nothing changes for the scaling and BvN decomposition algorithms with respect to the nonnegative case (Table 4 remains as it is for the general case).

Further investigations

Here, we give some further experiments to shed light into the behavior of the proposed preconditioners. We compare BvN * with BvN 𝑟 having the same number of permutation matrices. For this purpose, we created a BvN * preconditioner and used the number 𝑟 of its permutation matrices to create another preconditioner BvN 𝑟 by just taking the first 𝑟 permutation matrices from a BvN decomposition obtained by Algorithm 1. The results are shown in Table 7, where we give the sum of the coefficients used in constructing the preconditioners and the number of (F)GMRES iterations for a subset of matrices from Table 6, where BvN * led to a short running time.

As seen in Table 7, BvN 𝑟 has a larger sum of coefficients than BvN * with the same number 𝑟 of permutation matrices. This is expected as the Algorithm 1 finds coefficients in a non-increasing order; hence the first 𝑟 coefficients are the largest 𝑟 ones. The number of (F)GMRES iterations is usually inversely proportional to the sum of the coefficients; a higher sum of coefficients usually results in a smaller number of iterations for the same matrix.

Conclusions and open questions

We introduced a class of preconditioners for general sparse matrices based on the Birkhoff-von Neumann decomposition of doubly stochastic matrices. These preconditioners are aimed primarily at solving challenging linear systems with highly unstructured and indefinite coefficient matrices in parallel computing environments. We presented some theoretical results and numerical experiments on linear systems from a variety of applications. We regard this work to be a proof Table 5. The number of GMRES iterations with different number of permutation matrices in 𝑀 with tolerance 10 -6 . We ran GMRES (without restart) for at most 3000 iterations; the symbol "-" flags cases where GMRES did not converge to the required tolerance. Matrices have negative and positive entries.

BvN𝑟 it.

BvN* matrix ILU(0) Table 6. The complexity of the preconditioners and the running time of the solver. The complexities of BvN𝑟 are given as the ratio (nnz(𝐿 + 𝑈 ) -𝑛)/ nnz, where 𝐿 and 𝑈 are the factors of the preconditioner (ILU(0) has a complexity of 1.0), and that of BvN* is given as nnz(𝑀 )/ nnz(𝐴). The symbol "-" in the running time column flags the cases where (F)GMRES were unsuccessful. Matrices have negative and positive entries. of concept realization; many challenging questions remain to be investigated to render the proposed preconditioners competitive with the standard approaches. Based on our current theoretical findings, we suggest the use of proposed preconditioners within a Krylov subspace method. There are two ways to go about this. In the first one, the preconditioners are built to satisfy a sufficiency condition that we identified (Theorem 4). This way, the application of the preconditioner requires a small number of highly concurrent steps, where there is no data dependency within a step. In the second alternative, one obtains the LU decomposition of the preconditioners that are built arbitrarily. Here, the proposed preconditioner is therefore constructed as a complete factorization of an incomplete matrix. We demonstrated that using around eight matrices is good enough for this purpose. Beyond that number, the LU decomposition of the preconditioner can become a bottleneck (but remains always cheaper than that of the original matrix). Is there a special way to order these preconditioners for smaller fill-in?

The construction of the preconditioners needs an efficient doubly stochastic scaling algorithm. The known algorithms for this purpose are iterative schemes whose computational cores are sparse matrix-vector multiply operations whose efficient parallelization is well known. For constructing these preconditioners, we then need a bottleneck matching algorithm. The exact algorithms for this purpose are efficient on sequential execution environments, but hard to parallelize efficiently. There are efficiently parallel heuristics for matching problems [START_REF] Halappanavar | Codesign lessons learned from implementing graph matching on multithreaded architectures[END_REF], and more research is needed to parallelize the heuristic for obtaining a BvN decomposition while keeping an eye on the quality of the preconditioner.

Fig. 1 .

 1 Fig. 1. A greedy heuristic for obtaining a BvN decomposition.

  • • • 𝐴 2,𝑠 . . .

					⎤
		. . .	. . .	. . .	⎥ ⎥ ⎥ ⎦
	0	0	• • • 𝐴 𝑠,𝑠

Table 1 .

 1 Test matrices, their original size 𝑛 𝐴 and the size of the largest fully indecomposable block 𝑛. The experiments were done using the largest fully indecomposable block (𝑛 is the effective size in the following experiments). The last matrix's full name is FEM_3D_thermal1.

	matrix	𝑛 𝐴	nnz 𝐴	𝑛	nnz	𝑘
	slide	20191 1192535 19140 1191421	900
	two-dom	22200 1188152 20880 1186500	938
	watson5a	1854	8626	1765	6387	28
	watson4a	468	2459	364	1480	134
	bp_200	822	3802	40	125	47
	gemat11	4929	33108	4578	31425	234
	gemat12	4929	33044	4552	31184	260
	lns_3937	3937	25407	3558	24002	205
	mahindas	1258	7682	589	4744	722
	orani678	2529	90158	1830	47823 6150
	sherman2	1080	23094	870	19256	334
	west0655	655	2808	452	1966	307
	west0989	989	3518	720	2604	315
	west1505	1505	5414	1099	3988	345
	west2021	2021	7310	1500	5453	376
	circuit_3	12127	48137	7607	34024	322
	bayer09	3083	11767	1210	6001	300
	bayer10	13436	71594 10803	62238	4
	lhr01	1477	18427	1171	15914	473
	lhr02	2954	36875	1171	15914	476
	appu	14000 1853104 14000 1853104 1672
	raefsky4	19779 1316789 19779 1316789	774
	venkat25	62424 1717763 62424 1717763	375
	utm5940	5940	83842	5794	83148	291
	bundle1	10581	770811 10581	770811 9508
	fp	7548	834222	7548	834222 8008
	dw8192	8192	41746	8192	41746	155
	FEM_3D 17880	430740 17880	430740	329

  Performance profiles for the number of iterations with ILU(0), BvN𝑟 with 𝑟 = 2, 4, 8, and BvN*.

		0.5												0.5
		0.45												0.45
	matrix Fraction of test cases	0.1 0.15 0.2 0.25 0.3 0.35 0.4	ILU(0)	1		2			4		8 Fraction of test cases	0.1 0.15 0.2 0.25 0.3 0.35 0.4	16		32			64 𝑟	*	it.
	slide	0 0.05	1	1.1	-1102 1.2 1.3 ilu	1.4 BvN 2	1158 1.5 1.6 BvN 4	1.7	1.8 BvN 8	946 1.9 BvN *	2	787	0 0.05	567 1 1.1	1.2	566 1.3 1.4 ilu BvN 2	1.5	131 10 1.6 1.7 1.8 BvN 4 BvN 8	1.9 BvN *	2 884
	two-dom		1328	-	=	-			-		2774		2088		1719	=		1429 -	-
	watson5a (a) Nonnegative matrices -453 516			610		806 (b) General matrices 908 908			908 5	621
	watson4a		49 146		133			139		124		93		32			5 12	129
	bp_200 Fig. 2.			-	36		26			25		12		7		3			3 8	25
	gemat11		176	-		-		2857		2118		1406		664			169 -	-
	gemat12		263	-		-		2906		2237		1591		-			145 -	-
	lns_3937		349 2067		1286			666		229		56		14			5 10	1138
	mahindas			-303		283			245		173		87		40			21 7	278
	orani678			-371		360			334		303		283		264			221 5	335
	sherman2		19 324		238			157		96		51		21			7 9	235
	west0655		71 331		275			188		139		-		-			31 9	270
	west0989			-194		167			114		63		35		19			9 8	165
	west1505			-335		254			166		91		49		27			10 9	241
	west2021			-406		323			217		-		47		34			13 8	315
	circuit_3			-	-		-			-		-		-		-			1426 -	-
	bayer09		124 521		434			291		204		110		28			7 7	393
	bayer10			-	-		-			-		-		-		-			--	-
	lhr01				-574		350			216		130		86		42			18 7	313
	lhr02				-565		351			219		136		87		45			19 8	314
	appu			31	50		50			52		56		62		74			82 7	50
	raefsky4			-	-		2728		2735		2444		1559		1739			1230 6	2370
	venkat25			-	-		-			-		-		-		-			--	-
	utm5940			-2633		2115			-		-		-		153			34 7	2069
	bundle1		25 163		150			139		149		174		158			139 5	149
	fp				-	-		-			-		-		-		-			--	-
	dw8192			-	-		-		2398		1516		650		81			8 -	-
	FEM_3D			7	37		35			38		72		58		25			6 9	33
							Number of successful runs with all 90 matrices
				32	58		71			60		57		53		58			61 7	77

Table 4 .

 4 The total number of sparse matrix-vector and sparse matrix transpose-vector multiply operations (MVP) in the scaling algorithm, the running time of the scaling algorithm and the partial BvN decomposition algorithm in seconds.

			running time (.s)
	matrix	MVP Scaling	BvN
	two-dom	438	0.67	3.03
	watson4a	40368	0.44	0.00
	gemat11	2286	0.28	0.04
	gemat12	843216	1.01	0.04
	lns_3937	3732	0.27	0.04
	sherman2	1966	0.06	0.01
	west0655	1846	0.03	0.00
	bayer09	3626	0.10	0.01
	appu	46	0.15	1.57
	bundle1	116	0.13	0.78
	FEM_3D	52	0.04	0.81

Table 7 .

 7 Comparing BvN* with BvN𝑟 having the same number of permutation matrices in terms of the sum of the coefficient of permutation matrices and the number of (F)GMRES iterations (it.).

			Complexity of 𝑀				Running time
			BvN𝑟					BvN𝑟
		2	4	8	BvN* ILU(0)	2	4	8	BvN*
	slide	0.03	0.05	4.04	0.10 38.67	300.65	176.02	148.79	1017.34
	two-dom 0.03	0.04	1.31	0.08 56.40	183.60	-	205.53	1122.33
	watson4a 0.28	0.36	0.44	0.66	0.05	0.28	0.24	0.19	6.58
	gemat11 0.23	0.35	0.65	0.61 16.57	-	-	1548.86	-
	gemat12 0.21	0.31	0.58	0.61 34.15	-	-	1578.89	-
	lns_3937 0.29	0.55	1.66	0.63	4.68	186.25	-	3.12	359.80
	sherman2 0.06	0.10	0.14	0.26	0.02	0.80	0.25	0.09	18.03
	west0655 0.42	0.69	1.08	0.62	0.06	0.65	0.38	0.15	13.53
	bayer09	0.32	0.52	0.90	0.51	0.05	0.62	0.21	0.08	21.12
	appu	0.01	0.02	7.05	0.04	1.64	2.06	2.16	4.64	32.71
	raefsky4 0.02	0.02	0.29	0.07 206.02	262.48	116.39	182.10	1027.58
	venkat25 0.05	0.53	3.61	0.23 94.29	-	5462.64	-	20382.59
	bundle1	0.01	0.01	0.02	0.03	0.50	4.81	4.84	5.46	76.47
	FEM_3D 0.05	0.08	0.78	0.27	0.34	1.10	1.10	3.11	30.07

Matrices have negative and positive entries.

		∑	𝛼 𝑖	it.
	matrix	BvN*	BvN𝑟 BvN*	BvN𝑟
	watson4a	0.43	0.57	124	111
	lns_3937	0.29	0.80	881	68
	sherman2	0.24	0.61	226	59
	west0655	0.12	0.45	222	97
	bayer09	0.19	0.50	216	59
	appu	0.09	0.13	50	55
	bundle1	0.01	0.01	125	110
	FEM_3D	0.35	0.55	39	79

Acknowledgment:

The work of Michele Benzi was supported in part by NSF grant DMS-1418889. Bora Uçar was supported in part by French National Research Agency (ANR) project SOLHAR (ANR-13-MONU-0007). We thank Alex Pothen for his contributions to this work. This work resulted from the collaborative environment offered by the Dagstuhl Seminar 14461 on High-performance Graph Algorithms and Applications in Computational Science (November 9-14, 2014).