Yassine El Khadiri

Gabriel Corona

François Cédric Rose

Yassine El-Khadiri

François Charpillet

Sleep Activity Recognition using Binary Motion Sensors

Keywords: Inference mechanisms, Unsupervised learning, Senior citizens, Smart homes, Ambient Assisted Living

come   L'archive ouverte pluridisciplinaire

Introduction

The world population is converging towards 17% being over 60 years old by 2050 [START_REF] He | International population reports -an aging world: 2015[END_REF]. Care-taking of this senior population is getting more difficult given the general socio-economic context of the active and senior demographic. From the World bank database, the aged dependency ratio which is the ratio of people older than 64 to the working age population went from 8.6% in 1960 to 12.9% in 2016 [START_REF] Bank | Age dependency ratio, old[END_REF]. According to the United States Census Bureau, around 29% of adults aging 65 or older live alone [START_REF]The majority of children live with two parents, census bureau reports[END_REF]. The US national median cost of a semi-private Room in a nursing home is projected to grow by 243% from 2017 to 2047 [START_REF] Genworth | Cost of care survey 2017 -summary[END_REF]. This paper takes part in a project which goal is to develop an ambient health care monitoring system that minimizes the senior's involvement and interaction with it. The solution should then not need any type of learning or setup to start giving immediate accurate results. We present the algorithm we use for the Sleeping Activity of Daily Living.

Our method works on environments equipped with binary motion sensors. These sensors report motion that happens in their field of view. The sensors we use report a timestamp for when they detect movement and another one for when they no longer detect it. The reported activity data can then be represented by a boolean time-series. Using this representation has many advantages. It is easy to work with but also easy to translate into from more complex representations coming from non binary motion sensors like depth of field cameras for instance.

Related work

Formal sleep study in the medical field or Sleep medicine is done through Polysomnography [START_REF] Haba-Rubio | Evaluation instruments for sleep disorders: A brief history of polysomnography and sleep medicine[END_REF] which is a diagnostic tool in sleep medicine done in Sleep Labs. However, nowadays, there is a surge of Directto-consumer trackers and devices popularizing and bringing perhaps limited but still usefull sleep analysis tools to the masses [START_REF]A third of people track their health or fitness[END_REF][START_REF] Pogue | What fitbit's 6 billion nights of sleep data reveals about us[END_REF]. These tools come in different shapes and forms. Some are wearable sensors akin to smart watches [START_REF]Fitbit[END_REF][START_REF]Biostrap[END_REF]. Others are smartphone applications which come either as companions for the Wearables or as stand-alones [START_REF]Sleep cycle alarm clock[END_REF]. There is also more involved solutions like general public ambient smart-home sensors [START_REF]Sleeptracker[END_REF][START_REF]Beddit[END_REF] or whole remote Healthcare Monitoring systems [START_REF]Keystone care[END_REF][START_REF]Cormarch remote care for seniors[END_REF].

These applications make use of different models and methods to track and follow the user's sleep schedule.

Some are based on heartbeat monitoring [START_REF] Benedetto | Assessment of the fitbit charge 2 for monitoring heart rate[END_REF][START_REF] Coppetti | Accuracy of smartphone apps for heart rate measurement[END_REF] or accelerometers (actigraphybased solutions) [START_REF] Hyun-Ah | Comparison of wearable activity tracker with actigraphy for sleep evaluation and circadian rest-activity rhythm measurement in healthy young adults[END_REF], while others use simpler features like smartphone events (screen activations [START_REF] Abdullah | Towards circadian computing: "early to bed and early to rise" makes some of us unhealthy and sleep deprived[END_REF][START_REF] Cuttone | Sensiblesleep: A bayesian model for learning sleep patterns from smartphone events[END_REF], audio features [START_REF] Ren | Fine-grained sleep monitoring: Hearing your breathing with smartphones[END_REF], app usage [START_REF] Chen | Unobtrusive sleep monitoring using smartphones[END_REF]. . . ).

There are three different axes that need exploring to get an overview of these different methods for sleep monitoring: types of sensors, types of collected data, and types of algorithms.

These axes are not totally independent from one another, but they do offer a good framework for setting the context of any given work on this subject.

Following what kind of application is being sought for, sleep tracking can be done using different types of sensors. These can be categorized into two main groups: ambient or wearable. Perhaps the most wellknown ones are of the wearable kind with the big surge of smart-watches and activity tracker wristbands these past few years [START_REF]A third of people track their health or fitness[END_REF]. However, we are also noticing the same happening for the adoption of ambient smart-home solutions at least on the end consumer market [START_REF]Sleeptracker[END_REF]. Ambient sensors are in the middle ground between sleep medicine style studies using Sleep Labs and EEG setups which are aimed at technically advanced and specialized use cases and wearables which are mainly used outdoors during the day.

Ambient sensor setups offer a lot of possibilities in terms of data that can be collected. A simple and straight forward one would be video capture which is often frowned upon and not desirable for obvious privacy concerns. In the context of home automation, motion sensors (PIR sensors) are the most used data sources for activity tracking indoors [START_REF] Luo | Simultaneous indoor tracking and activity recognition using pyroelectric infrared sensors[END_REF]. These can also come packaged with temperature, light, humidity and other types of sensors in general [START_REF]Fibaro motion sensors[END_REF]. Other, perhaps more exotic, sensors allow for capturing 3d body skeletons and silhouettes [START_REF] Dubois | A gait analysis method based on a depth camera for fall prevention[END_REF] or heat signatures with radio signals [START_REF] Hsu | Zero-effort in-home sleep and insomnia monitoring using radio signals[END_REF].

Algorithms for sleep or ADL tracking in general come in two groups [START_REF] Rashidi | A survey on ambient-assisted living tools for older adults[END_REF]: supervised and unsupervised. Some well-known examples of supervised methods used for human activity recognition are Classifiers like Naive Bayes and Random Forest [START_REF] Nef | Evaluation of three state-of-the-art classifiers for recognition of activities of daily living from smart home ambient data[END_REF], Hidden Markov Models [START_REF] Lu | Activity recognition in smart homes[END_REF] and Deep Convolutional Neural Networks [START_REF] Yang | Deep convolutional neural networks on multichannel time series for human activity recognition[END_REF]. Even though the majority of activity recognition or more specifically ambient activity recognition algorithms are supervised, they are inherently hard to scale given not only to their need for annotated training data but also due to inherent differences in the lifestyles of people. Unsupervised methods are mostly model based inference algorithms that can work on a small amount of data. They range from simple rule-based algorithms [START_REF] Abdullah | Towards circadian computing: "early to bed and early to rise" makes some of us unhealthy and sleep deprived[END_REF] to Monte-Carlo Markov Chain setups [START_REF] Cuttone | Sensiblesleep: A bayesian model for learning sleep patterns from smartphone events[END_REF].

We propose an ambient assisted living approach that minimizes the patient's involvement with the monitoring system. We make use of motion sensors given their simplicity and availability in conjunction with a model based inference algorithm.

Method

Dataset

We use one of the many publicly available datasets of the WSU CASAS Project [START_REF] Cook | Learning setting-generalized activity models for smart spaces[END_REF] for the experimentation and evaluation of our algorithm. The CASAS project's purpose is to serve as a platform for the collection of real world data for research projects related to smart environments and health monitoring and intervention. The one we used in particular contains data collected from the home of a volunteer adult. The home is equipped with motion sensors that report binary ON and OFF data with a timestamp at each time they sense and stop sensing activity.

which means that we basically have access to a boolean time-series of sensor activations.

The dataset we use has a total of 219 nights worth of data spanning from 2010-11-04 to 2011-06-11 with over 1,713,128 sensor activations reported from the sensors placed all around the resident's house which floor-plan.

Figure 1 shows reported data from these sensors centered around midnight and re-sampled on slices of 10 minutes. As we can see there is an obvious change of rate in sensor activations when the resident is asleep at night. This is the main assumption we make for the elaboration of our algorithm for finding the sleeping periods which we will assume is the interval where we have different sensor activations (in most cases fewer).

The input data for the algorithm for each night inference is centered around midnight and re-sampled on slices of 10 minutes.

Algorithm

Model

The main assumption being that the reported sensor data is strongly correlated with the level of activity in the apartment. Periods with very few sensor activations would correlate to calmer activities like Reading or Relaxing, and even longer periods to Sleeping. We also assume that the resident is awake at noon which is most likely the case especially for elder senior residents in nursing homes for example.

Our method is based on the Bayesian Switch-point Model [START_REF] Maguire | The time intervals between industrial accidents[END_REF], sometimes also known by the name of Bayesian Change Point Analysis: given a sample of observations from a process of randomly occurring events over a period of time, we search for a change in the rate at which events occur. For example an increase or decrease in the frequency of the observations.

We name the period of low activity the Sleep period, and the period of high activity the Wake period.

Since sensor observations (M t ) are boolean time series, they are modeled after a Bernoulli law (1) of parameter P t : the probability that the sensor would activate at the slice t. This parameter changes according the period the slice t is in: Let P s be the probability of observing a sensor activation during the Sleep period and P w of observing a sensor activation during the wake period. Since the goal is to find when the Sleep period starts and ends we introduce the two variables T s and T w that correspond to its boundaries. This enables us to build a boolean variable S t (2) that conveniently reports which period each time slice t is in.

M t ∼ B(P t ), P t = P s if S t = 0 P w if S t = 1 (1)
S t = 0 if T s ≤ t < T w 1 if t < T s or t ≥ T w (2)

Maximum a posteriori estimate

We consider that we have a stronger prior knowledge of the wake times T w and sleep duration D = T w -T s than on the sleep time T s . This way the value of T s is not too much constrained towards its prior.

Inferring T s and T w is done by computing the maximum a posteriori (MAP) estimate P of T s = t s , T w = t w , P s = p s , P w = p w with [t s , t w ] being the sleeping period where a sensor get activated with probability p s and [0, t s [∪]t w , N ] the wake period where a sensor gets activated with probability p w given the observations of sensor activations M = {M t , ∀t ∈ [0, N ]}:

T s , T w = argmax
P (M | T s , T w , P s , P w ) being the likelihood of T s = t s , T w = t w , P s = p s , P w = p w given the observations M .

We consider that P s and P w have a uniform prior meaning that P (P s , P w ) cancels out in the maximization process and we have:

P (T s = t s , T w = t w ) = d P (T s , T w , D = d) =P (T s , T w , D = t w -t s ) + d =tw-ts P (T s , T w , D) =P (T w = t w ) × P (D = t w -t s ) (5)
As for the T w and D priors we choose a normal distributions centered around 8:00 AM in the morning for T w and 8 hours of sleep for the night duration and both with a standard deviation of an hour and 30 minutes:

T w ∼ N (8 : 00 AM, 1h30) D ∼ N (8 hours, 1h30) (6)
In practice the T w and D values are discretized. Figure 2 is a graphical representation of the model described above. The slate t = 0, . . . , N -1 represents the number of time slices that constitutes the data record spanning from noon of a day d to noon of day d + 1. The slate R ∈ Rooms represents the number of sensor timeseries available.

From equation 3 and 4 we have to maximize the posterior:

P (T w )P (D) prior t P (M t | T s , T w , P s , P w ) likelihood L (7)
In order to maximize the posterior we have to maximize L over P s and P w : max

Soft evidence

However we have to take into consideration that some sensor observations in some rooms don't give any information related to our model parameters. For example observing a lack of movement doesn't systematically correlate with a passive activity anymore since the resident could be active elsewhere. In other words, our evidence on the observations is no longer certain when taking into consideration all the different sensors scattered in the different rooms. Thus we need to introduce soft-evidence in the inference calculations.

Soft-evidence can be modeled over the observations M by the following distribution: e(M t | X t ) [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF]. It isn't a model parameter but an external process we hook into the MAP computation to model the uncertainty of the observations. So we get to freely set the values of e:

e M t = 0 M t = 1 X t = 0 0.3 0.7 X t = 1 0 1
The new objective function will be the expectation of the log-posterior:

M ∼e(X) log (P (T w ) × P (D) × L) de(X) (12) i.e.
log P (T w ) + log P (D) + log L

For a given T s and T w , the optimum value of P s is the one that maximizes: log P (M t|St=0 | T s , T w , P s , P w )de(X)

= log P t E[M t ] + log(1 -P ) t (1 -E[M t ]) Which is: Ps = t|St=0 E[M t ] Card({t|S t = 0}) (14)
The same goes for P w .

Method evaluation

We evaluate our algorithm using the sleep annotations present in the dataset. We run the inference for each night then compute the error over the inferred values of T s and T w : We also evaluate it as a binary classification problem for each time slice by constructing the vectors of annotated times S = (S t | t ∈ 0..N ) and inferred times Ŝ = ( Ŝt | t ∈ 0..N ). We can thus compute the accuracy, precision, recall and F1 scores over all the nights.

e x = T x -T annotated x

Results

The results on the annotated dataset yields a global accuracy 0.94 with an F1 score of 0.90. The precision and recall scores are respectively 0.85 and 0.95. Figure 4 shows a more detailed view for these results on an histogram over all of the dataset's nights. From the accuracy histogram, 80.82% of the evaluated nights have an accuracy over 0.9. 44.74% for the precision, 89.5% for the recall and 56.16% for the F1 score.

Discussion

The main contribution of this work is a method for sleep inference which is computationally inexpensive and applicable to all kinds of ambient sensor setups as long as the reported data holds motor activity information. This method is already directly usable on any new or existing installation as is or with minimum tuning on the T s and D priors. The use of soft evidence makes it possible to work on data without having to take into consideration the resident's whereabouts and possibly have in the mix sensor data not following the sleep-wake model without it impacting too much the inference results.

The accuracy, precision, recall, and F1 global and histogram scores show that this method is also in line with and even slightly better than the expected results from an unsupervised one [START_REF] Cuttone | Sensiblesleep: A bayesian model for learning sleep patterns from smartphone events[END_REF].

SensibleSleep [START_REF] Cuttone | Sensiblesleep: A bayesian model for learning sleep patterns from smartphone events[END_REF] is also a discrete event based method for extracting sleep patterns. The paper uses this method in conjunction with smartphone events. The inference is based on the Markov Chain Monte Carlo method [START_REF] Andrieu | An introduction to mcmc for machine learning[END_REF] that fits the model to the observed data and outputs the posterior distribution of the parameters. In contrast, we show that by using our model, directly computing a closed form solution yields comparable results while being computationally inexpensive even on smaller bin sizes (our method discretizes the day into 10 minute bins in-stead of 15 minutes for theirs).

Another advantage of this method is the fact that we do not assume P s ≤ P w in the inference algorithm since it is not necessarily true all the time. For example chances are that we observe more movements in the bedroom during the night than the day either because of sleep agitation or bathroom visits or sensitive sensors. However this explains outlier cases where the algorithm finds a segment of higher likelihood with more sensor activations than the true sleep period.

There are also other limits to how reliable this algorithm can be, since as shown on figure 3, the inference seems inclined into overestimating wake times and underestimating sleep times.

A possible enhancement to the algorithm can be to add a learning component on the priors so the more data gets collected the more they get refined. Another one might be a way to compute the confidence in the inferred sleep periods as a starter for work on alerts of sleep agitation or deprivation or unusual and abnormal activity.

However this algorithm in its current state is already deployed in 16 installations including nursing homes and private apartments. The wake time inference is used by the nursing staff as a way to approximate the time when residents wake up so they can go help them in their morning routine. From the inferred sleep period we can also infer nocturnal activity for potentially fragmented sleep problems or medical problems related to frequent bathroom visits.

Conclusion

We have presented an unsupervised method for sleep period detection using binary time-series of motor activity data. This method is quite simple and computationally inexpensive which makes it possible to integrate into existing medical monitoring devices. It is also cost-effective as it doesn't need any special setup (calibration, learning phase or special sensor placement) that for example depends on the geometry or the nature of the residency's different rooms. The method also supports a dynamically changing number of sensors. This algorithm is currently being

Figure 1 :

 1 Figure 1: Illustration of sensor data centered around midnight showing a clear change in the rate of sensor activations.

  s , T w , P s , P w | M ) (3) with: P (T s , T w , P s , P w | M ) = P (M | T s , T w , P s , P w ) × P (T w ) × P (D)

Figure 3 :

 3 Figure 3: Boxplot of the error on T s , T w , and D

Figure 3

 3 Figure 3 contains the boxplots of the distributions of the relative error on T s , T w , and the night duration D rounded on hours. Both medians are centered around 0 minutes. e s whiskers span from -3 hours and 10 minutes to +2 hours while e w 's span from -1 hour 50 minutes to +2 hours 50 minutes.

Figure 4 :

 4 Figure 4: Histogram of the computed accuracy, precision, recall and F1 scores for each night

used in nursing home rooms and apartments for the experimentation phase of the project this paper has originated from.


Emmanuel Agullo
email: emmanuel.agullo@inria.fr

Eric Darve
email: eric.darve@stanford.edu

AND Yuval Harness
email: yuval.harness@inria.fr

Low-rank Factorizations in Data Sparse Hierarchical Algorithms for Preconditioning Symmetric Positive Definite Matrices

Keywords: Preconditioning, Symmetric Positive Definite, Data Sparse, Hierarchical Algorithms, Low-rank Factorization, Minimal Condition Number AMS subject classifications. 15A16, 15B99, 65F08, 65F30, 65F35, 65F50

We consider the problem of choosing low-rank factorizations in data sparse matrix approximations for preconditioning large scale symmetric positive definite matrices. These approximations are memory efficient schemes that rely on hierarchical matrix partitioning and compression of certain sub-blocks of the matrix. Typically, these matrix approximations can be constructed very fast, and their matrix product can be applied rapidly as well. The common practice is to express the compressed sub-blocks by low-rank factorizations, and the main contribution of this work is the numerical and spectral analysis of SPD preconditioning schemes represented by 2 × 2 block matrices, whose off-diagonal sub-blocks are low-rank approximations of the original matrix off-diagonal subblocks. We propose an optimal choice of low-rank approximations which minimizes the condition number of the preconditioned system, and demonstrate that the analysis can be applied to the class of hierarchically off-diagonal low-rank matrix approximations. Spectral estimates that take into account the error propagation through levels of the hierarchy which quantify the impact of the choice of low-rank compression on the global condition number are provided. The numerical results indicate that the properties of the preconditioning scheme using proper low-rank compression are superior to employing standard choices for low-rank compression. A major goal of this work is to provide an insight into how proper reweighted prior to low-rank compression influences the condition number for a simple case, which would lead to an extended analysis for more general and more efficient hierarchical matrix approximation techniques.

Introduction.

In this paper we consider preconditioning for iterative solution of large scale linear systems [START_REF] Agullo | Nearly optimal fast preconditioning of symmetric positive definite matrices[END_REF] Ax = b , where A ∈ R n×n is a symmetric positive definite (SPD) matrix. Such systems arise in a wide range of engineering applications, as means to model and understand physical phenomena. Typical example is the result of a finite element discretization of underlying differential equations of a boundary value problem. In many practical applications the matrix A becomes ill-conditioned and, thus, challenging for iterative techniques. In that case the use of preconditioned iterative methods, such as the preconditioned conjugate gradient (PCG) [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF][START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] technique, becomes an imperative.

The choice of a suitable preconditioning scheme can, often, drastically improve the convergence behavior of the iterative method and, generally, plays a vital role in the success of solving the system.

A preconditioning scheme for linear systems is, essentially, composed of linear operations or matrices that approximate A -1 (1), but with considerable less computational effort than explicitly inverting A. Transforming the system (1) with such a scheme is called the preconditioned system. The major concern when setting up a preconditioning scheme is to ensure that the preconditioned system has a bounded condition number, and that the number of iterations to convergence in an iterative method remains small while maintaining low associated complexity and reduced mem-ory cost. The literature on preconditioning techniques is vast, and many robust and efficient methods have been introduced in the last 50 years. These include, among others, incomplete factorization schemes such as ILU and Incomplete Cholesky, sparse matrix approximations, polynomial techniques, domain decomposition methods, as well as multigrid and algebraic multilevel iterations schemes. For a recent comprehensive review on this topic see [START_REF] Wathen | Preconditioning[END_REF].

The main contribution of this work is the numerical and spectral analysis of SPD preconditioning schemes represented by 2 × 2 block matrices, whose off-diagonal subblocks are low-rank approximations of the original matrix off-diagonal sub-blocks. We re-examine the way low-rank factorizations are obtained, by considering reweighting of the sub-blocks prior to the low-rank compression. Reweighting can be done in many ways, e.g., diagonal scaling, and the fundamental question that we attempt to answer is: which reweighting is optimal with respect to the condition number of the preconditioned system?

The mathematical theory for 2 × 2 matrices is derived in section 2. We present an optimal 1-level preconditioning scheme using proper reweighting prior to low-rank compression, which minimizes the spectral condition number of the preconditioned system. Thus, a preconditioning scheme employing such low-rank factorizations is expected to attain the same condition number with less computational resources and associated complexity, compared to employing other standard techniques for the lowrank factorizations. Spectral analysis shows that the scheme maps both small and large eigenvalues of the original system exactly to 1. This feature is of great importance to Krylov subspace methods, since it is equivalent to the minimization of the effective degree of the minimal polynomial of A that defines the maximal dimension of the search space.

In section 3 we propose an application of the 1-level theory for hierarchically off-diagonal low-rank (HODLR) matrix structure, as means to demonstrate the applicability of the 1-level theory to the hierarchical multilevel case. We also provide spectral estimates that take into account the error propagation through levels of the hierarchy. This leads to quantification of the impact of the reweighting on the global condition number of the preconditioned system. In essence, weighted HODLR locally minimizes the condition number at each level of the hierarchy by approximately filtering the smallest and largest eigenvalues. Since this approach is employed hierarchically, it effectively creates a strong effect of global spectrum clustering.

The HODLR structure is a member of a wide class of hierarchical data sparse approximations. These approximations rely on the fact that the matrix can be subdivided into a hierarchy of smaller block matrices, and certain sub-blocks can be efficiently approximated as low-rank matrices by low-rank factorizations. The low-rank compressions of sufficient sub-blocks leads to a dramatic reduction of the complexity and computational cost. The best known example for such schemes is the class of hierarchical matrix (H-matrix) approximations [START_REF] Grasedyck | Construction and Arithmetics of H-Matrices[END_REF][START_REF] Hackbusch | A Sparse Matrix Arithmetic Based on H-Matrices. Part I: Introduction to H-Matrices[END_REF][START_REF] Hackbusch | A Sparse H-Matrix Arithmetic. Part II: Application to Multi-Dimensional Problems[END_REF][START_REF] Bebendorf | Hierarchical Matrices: A Means to Efficiently Solve Elliptic Boundary Value Problems[END_REF] which has gained growing attention in recent years.

To the contrary of the more general strong hierarchical matrix structure which allows further decomposition of the off-diagonal blocks into low-rank and full-rank blocks, HODLR is a weak hierarchical matrix structure, which relies on a single lowrank compression for the off-diagonal blocks. Closely related to HODLR is the hierarchically semi-separable (HSS) [START_REF] Chandrasekaran | A Fast U LV Decomposition Solver for Hierarchically Semiseparable Representations[END_REF][START_REF] Xia | Fast algorithms for hierarchically semisep-arable matrices[END_REF] structure , which is, in fact, a HODLR matrix format possessing a nested off-diagonal low-rank structure.

Essentially, weak hierarchical methods, i.e., HODLR and HSS, are not considered competitive with the more general strong hierarchical matrix methods, when the un-This manuscript is for review purposes only.

derlying problem is of very large scale. However, the proposed study provides a novel theoretical basis for optimality conditions of hierarchical preconditioning schemes.

Thus, the presented analysis can serve as starting point for a more general theory on optimal H-matrix preconditioning which is deferred to future work.

The weighted HODLR scheme is similar in nature to the methods proposed in [START_REF] Xia | Robust Approximate Cholesky Factorization of Rank-Structured Symmetric Positive Definite Matrices[END_REF][START_REF] Xia | Effective and Robust Preconditioning of General SPD Matrices via Structured Incomplete Factorization[END_REF] which propose practical HSS schemes that rely on similar ideas of reweighting prior to compression, but without the complete numerical and spectral analysis of this study. The costs to apply these multilevel preconditioners are about O(n), where n is the matrix size.

The experimental part of this work, whose goal is to demonstrate the effectiveness of properly chosen reweighting for the low-rank approximations, is given in section 4.

The section contains a detailed comparative study of HODLR preconditioning using different methods for the low-rank compressions. As alternatives to the proper reweighting strategy, we consider the conventional low-rank approximation in the 2-norm and the low-rank approximation with constraints [START_REF] Bebendorf | Hierarchical matrix approximation with blockwise constraints[END_REF]. The latter employs lowrank approximations that also preserve constraints, forcing sub-blocks of the preconditioning scheme to be identical to the corresponding sub-blocks of the input matrix on predetermined subspaces. Employing the method for preconditioning SPD matrices of discretized elliptic PDEs has been demonstrated in [START_REF] Bebendorf | On the spectral equivalence of hierarchical matrix preconditioners for elliptic problems[END_REF], and a similar approach for non-symmetric sparse matrices has been recently suggested in [START_REF] Yang | Sparse Hierarchical Solvers with Guaranteed Convergence[END_REF].

The numerical results indicate, that employing proper reweighting prior to lowrank compression, leads to a HODLR preconditioning scheme that requires far less computational resources for the same quality of convergence performance compared to using other low-rank compression techniques. The experiments also show, that the HODLR preconditioning scheme with proper reweighting retains the SPD property of the system when other standard techniques fail, and remains efficient and robust even if low accuracy compression is employed with ranks of O(1) for the low-rank approximations of the sub-blocks. Summary and plans for future work follow in section 5.

2.

The Optimal One-level Preconditioning Scheme. In this section we introduce the optimal 1-level scheme for the preconditioning of SPD matrices. We consider an input n×n SPD matrix A with a 2×2 block structure and a corresponding 1-level approximation K,

(2)

A = A 1 M M T A 2 , K = A 1 U 1 SV T 2 V 2 SU T 1 A 2 , A i ∈ R ni×ni ,
where n = n 1 + n 2 , and the off-diagonal blocks of K are low-rank factorizations satisfying (3)

U 1 ∈ R n1×r , S ∈ R r×r , V 2 ∈ R n2×r ,
with a, typically, small rank r. The matrix U 1 is the interpolation operator, the matrix V 2 is the anterpolation operator, and the matrix S whose rank is r is known as the interaction operator. In some cases S is omitted, i.e., equivalently represented by an r × r identity matrix.

We present an explicit formula for a 1-level approximation, K (2), which minimizes the spectral condition number

(4) cond 2 R -T AR -1 = R -T AR -1 2 • RA -1 R T 2 ,
This manuscript is for review purposes only.

of the preconditioned system, (5)

R -T AR -1 x = R -T y ,
for any given rank r = 0, 1, . . . , min{n 1 , n 2 }, where • 2 is the 2-norm, and R denotes any square root (not necessarily principal) of K in the sense that ( 6)

K = R T R ∈ R n×n .
The key idea is to reweight the off-diagonal blocks prior the low-rank factorization.

Proper choice of reweighting leads to a minimum spectral condition number of the preconditioned system as well as clustering of the spectrum of the preconditioned system around 1.

We begin in subsection 2.1 by introducing the method for obtaining the minimum condition number low-rank approximation. In subsection 2.2 we provide the theorem on the minimum condition number property, including a detailed description of the spectral properties of the preconditioned system. A rigorous and detailed proof of the theorem is given in subsection 2.3.

2.1. Explicit Formula of the Optimal One-level Scheme. The construction of the minimum condition number 1-level preconditioner K (2) subject to

(7) rk U 1 SV T 2 ≤ r ,
is based on the following two-step method ensuring that the preconditioned matrix R -T AR -1 also inherits the SPD property of A:

1. Apply a two-sided block Jacobi transformation,

(8) R -T 1 0 0 R -T 2 • A • R -1 1 0 0 R -1 2 = I 1 R -T 1 M R -1 2 R -T 2 M T R -1 1 I 2 ,
where I i denotes the n i × n i identity matrix, and R i ∈ R ni×ni denotes a square root of

A i i.e., R T i R i = A i .
2. Extract the off-diagonal triple products (3) by setting,

U 1 = R T 1 U r , S = Σ r = diag(σ 1 , . . . , σ r ) , V 2 = R T 2 V r , (9)
where U r and V r are composed of the first r left and right, respectively, singular vectors of the singular value decomposition (SVD),

R -T 1 M R -1 2 = UΣV T , Σ = diag(σ 1 , . . . , σ min{n1,n2} ) .
The theory presented in this study implies, that in practice for a given rank bound r ≥ 0, any low-rank factorization, U 1 SV T 2 , satisfying

R -T 2 M T R -1 1 -R -T 2 U 1 SV T 2 R -1 1 2 ≤ R -T 2 M T R -1 1 -U r Σ r V T r 2 ,
would achieve the same minimal spectral condition number. However, the truncated SVD of the reweighted off-diagonal block, U r Σ r V T r , also ensures that the spectrum of the preconditioned system is optimally clustered. This observation is discussed and explained in the next subsection.

This manuscript is for review purposes only.

Minimal Condition Number and Spectral Analysis. Let us now focus

on the spectral properties of the preconditioned system, R -T AR -1 , where R is a square root of K (2) whose off-diagonal low-rank blocks are given by [START_REF] Chevalier | PT-Scotch: A tool for efficient parallel graph ordering[END_REF]. First, let us consider the degenerate case r = 0. In this case U 1 SV T 2 = 0 and the square root of K reduces to the following block diagonal form, ( 10)

R(r = 0) = R 1 0 0 R 2 .
The preconditioning scheme (5) with r = 0 is, in fact, the two-sided block Jacobi [START_REF] Chandrasekaran | A Fast U LV Decomposition Solver for Hierarchically Semiseparable Representations[END_REF].

There is a known result [START_REF] Elsner | A note on optimal block-scaling of matrices[END_REF] showing that the two-sided block Jacobi preconditioner ( 8) is optimal, in the sense that

cond 2 R -T 1 0 0 R -T 2 A R -1 1 0 0 R -1 2 ≤ cond 2 B T 1 0 0 B T 2 A B 1 0 0 B 2 ,
for any non-singular B1 0 0 B2

with the same dimensions and partition as R(r = 0) [START_REF] Davis | The University of Florida Sparse Matrix Collection[END_REF]. The analysis we present, thus, naturally extends this classic result.

The main results for the general case r ≥ 0 are presented in Theorem 1, whose principal component is the spectral analysis of the preconditioned system. Our proof shows that the spectrum of the two-sided block Jacobi preconditioned system [START_REF] Chandrasekaran | A Fast U LV Decomposition Solver for Hierarchically Semiseparable Representations[END_REF] contains (or equals to)

1 + σ 1 , . . . , 1 + σ min{n1,n2} , 1 -σ min{n1,n2} , . . . , 1 -σ 1 ,
where 1 -σ 1 and 1 + σ 1 are the smallest and largest, respectively, eigenvalues of the preconditioned system. Thus, the two-sided block Jacobi redistributes the spectrum of the matrix symmetrically around 1. We show that the optimal 1-level preconditioning scheme does the same, but also maps the largest r eigenvalues (1 + σ 1 , . . . , 1 + σ r ) and the smallest r eigenvalues (1 -σ r , . . . , 1 -σ 1 ) of ( 8) exactly to 1. Hence, the spectral condition number (4) as a function of r is

cond 2 R -T AR -1 = 1 + σ r+1 1 -σ r+1 .
An illustration of the spectral clustering done by the optimal 1-level preconditioning scheme is displayed in Figure 1.

Theorem 1. Let A = A 1 M M T A 2 , K = A 1 U 1 SV T 2 V 2 SU T 1 A 2 ,
have the same dimensions and partition where A is SPD, and let R i denote a square

root of A i , i.e., A i = R T i R i .
If the off-diagonal triple product approximation U 1 SV T 2 satisfies

U 1 = R T 1 U r , S = Σ r , V 2 = R T 2 V r ,
where U r and V r are composed of the first r left and right, respectively, singular vectors of the SVD,

(11) R -T 1 M R -1 2 = UΣV T , Σ = diag(σ 1 , . . . , σ min{n1,n2} ) , then:
This manuscript is for review purposes only. The spectrum of some SPD matrix A and the transformations it goes after preconditioning by block Jacobi (BJ) and the optimal 1-level preconditioning scheme are displayed. The spectra are ordered from the left to the right starting from A, followed by BJ and end up with the optimal scheme.

A BJ Weighted

1. The matrix K is SPD and possesses a square root, R, i.e., K = R T R.

2.

For any r < min{n 1 , n 2 } and any square root R, the spectrum of the preconditioned system is contained in ]0, 2[ and equal to

1 + σ r+1 , . . . , 1 + σ min{n1,n2} , 1, 1 -σ min{n1,n2} , . . . , 1 -σ r+1 .
3. Any inverse square root of K, R -1 , is a minimizer of the spectral condition number (4) in the sense that

cond 2 R -T AR -1 ≤ cond 2 R -T A R -1 , R T R = K ,
for any SPD matrix with the same dimensions and partition as K of the form,

K = A 1 M M T A 2 ,
whose off-diagonal blocks satisfy rk( M ) ≤ r.

Proof of Theorem 1.

The proof of Theorem 1 regarding the spectral properties relies on Lemma 1, while the minimum condition number property is based on the Cauchy (eigenvalue) interlacing theorem [23, p. 202]. The latter asserts that the eigenvalues of any principal submatrix of a symmetric matrix interlace those of the symmetric matrix. To be precise, if H ∈ R n×n is a partitioned symmetric matrix of the following form

H = E F F T G ,
This manuscript is for review purposes only.

in which E is a r × r principal submatrix, then for each i = 1, . . . , r,

λ i (H) ≥ λ i (E) ≥ λ i+n-r (H) , (12)
where the eigenvalues of the symmetric matrix H are assumed to be arranged in a decreasing order: ) where I i is the n i × n i identity matrix and δ ∈ R, and let σ 1 , . . . , σ min{n1,n2} denote the singular values of M.

λ 1 (H) ≥ λ 2 (H) ≥ • • • ≥ λ n (H) . Lemma 1. Let H = δI1 M M T δI2 ∈ R (n1+n2)×(n1+n2
1. If n 1 = n 2 then spec(H) = {δ -σ 1 , . . . , δ -σ n1 , δ + σ n1 , . . . , δ + σ 1 } . 2. If n 1 = n 2 then spec(H) = {δ -σ 1 , . . . , δ -σ min{n1,n2} , δ + σ min{n1,n2} , . . . , δ + σ 1 } ∪ {δ} ,
where the multiplicity of the eigenvalue δ is at least

|n 1 -n 2 |.
Proof. of Lemma 1.

Let us assume without the loss of generality that n 1 ≥ n 2 = m and let

M = UΣV T , U ∈ R n1×n2 , V ∈ R n2×n2 , denote the SVD of M. Let (13) W = 1 √ 2 U U V -V ∈ R n×n ,
whose blocks are given by

U = U if n 1 = n 2 U √ 2U ⊥ if n 1 > n 2 , V = V if n 1 = n 2 V 0 if n 1 > n 2 ,
where U ⊥ is an n 1 × (n 1 -n 2 ) matrix with orthonormal columns, whose range is orthogonal to the range of U,

U T U ⊥ = 0 ∈ R n2×(n1-n2) .
Direct calculations show that W is an orthonormal matrix satisfying

W T 0 M M T 0 W = S 1 0 0 -S 2 ,
where S i = diag σ 1 , . . . , σ min{n1,n2} , 0, . . . , 0 ∈ R ni×ni . Thus, by the orthogonality of W we obtain

W T δI 1 M M T δI 2 W = δI 1 + S 1 0 0 δI 2 -S 2 .
Hence, the spectrum of H is given by spec

(H) = {δ -σ 1 , . . . , δ -σ min{n1,n2} , δ + σ min{n1,n2} , . . . , δ + σ 1 } ∪ {δ} ,
where the multiplicity of δ is at least n 1 -n 2 . Note that in case n 2 > n 1 , we can simply interchange the principal blocks of H by reordering the columns and rows of H, and repeat the proof.

This manuscript is for review purposes only.

Proof. of Theorem 1.

Let K be a partitioned SPD matrix with the same dimensions and partition as K (2)

whose off-diagonal blocks rank is bounded by r,

K = A 1 U 1 SV T 2 V 2 SU T 1 A 2 , rk(U 1 SV T 2 ) ≤ r . If (λ, ζ) ∈ R × R n is an eigenpair of the preconditioned matrix (14) R -T A R -1 , K = R T R ,
then by employing the change of variables, ζ = Rξ, we obtain

R -T A R -1 ζ = λζ ⇔ R -T Aξ = λ Rξ ⇔ R -1 R -T Aξ = λξ . Since R -1 R -T = K -1
, we conclude that regardless to the particular choice of square root, R, the spectrum of the preconditioned system ( 14) remains unchanged.

Let R i ∈ R ni×ni denote a, generally, non-symmetric square root of A i . By direct calculations we obtain

(15) R -T 1 0 0 R -T 2 A R -1 1 0 0 R -1 2 = I 1 R -T 1 M R -1 2 R -T 2 M T R -1 1 I 2 ,
and by Lemma 1, the spectrum of ( 15) is contained in (or equal to)

{1 + σ 1 , . . . , 1 + σ min{n1,n2} , 1, 1 -σ min{n1,n2} , . . . , 1 -σ 1 } , where σ i are the singular values of R -T 1 M R -1 2 .
Since R i are non-singular, the preconditioned matrix ( 15) is SPD. Hence, we have

1 -σ 1 > 0 ⇒ spec I 1 R -T 1 M R -1 2 R -T 2 M T R -1 1 I 2 ⊂]0, 2[ .
Consider the specific choice of inverse square root of K,

R -1 = R -1 1 0 0 R -1 2 W D -1/2 W T , D = I 1 + S 1,r 0 0 I 2 -S 2,r .
The matrix W is orthogonal of the same form as [START_REF] Frieze | Fast Monte-carlo Algorithms for Finding Low-rank Approximations[END_REF] in the proof of Lemma 1 with respect to [START_REF] Grasedyck | Construction and Arithmetics of H-Matrices[END_REF], and S i,r = diag (σ 1 , . . . , σ r , 0, . . . , 0) ∈ R ni×ni where σ i are the

singular values of R -T 1 U 1 SV T 2 R -1 2 .
Setting the choice R -1 into (4) and employing the fact that the 2-norm is invariant under unitary transformations, we obtain

cond 2 ( R -T A R -1 ) = D -1/2 H D -1/2 2 D 1/2 H -1 D 1/2 2 ,
where H is an SPD matrix given by

H = W T WDW T W , D = I 1 + S 1,m 0 0 I 2 -S 2,m .
The matrices W (13) and S i,m are defined and constructed in the proof of Lemma 1.

Note that like W, the matrix W is orthogonal. Hence, the product W T W is also an orthogonal matrix.

Our definitions so far indicate that the following diagonal matrices,

D = I 1 0 0 (I 2 -S 2,r ) , D = (I 1 + S 1,r ) 0 0 I 2 . bound the diagonal matrix D D ≤ D ≤ D
in the sense that D -D and D -D are non-negative definite. Thus, applying the change of variables ξ = D -1/2 x and exploiting the properties of the Rayleigh quotient, we can write

D -1/2 H D -1/2 2 = max x =0 x T D -1/2 H D -1/2 x x T x = max ξ =0 ξ T Hξ ξ T Dξ ≥ max ξ =0 ξ T Hξ ξ T Dξ = max y =0 y T D -1/2 HD -1/2 y y T y = D -1/2 HD -1/2 2 , ( 16
)
where y = D 1/2 ξ. Using the same arguments it can also be shown that

D 1/2 H -1 D 1/2 2 ≥ D 1/2 H -1 D 1/2 2 .
Let Z = span{e r+1 , . . . , e n } where e i denotes the i-th canonical basis vector, and let P Z denote the orthogonal projection matrix on Z. The structure of D implies

D -1/2 HD -1/2 2 ≥ max P Z y =0 (P Z y) T D -1/2 HD -1/2 P Z y (P Z y) T P Z y = max P Z y =0
(P Z y) T HP Z y (P Z y) T P Z y .

Essentially, P T Z HP Z represents an (n -r) × (n -r) principal block of a unitarily equivalent matrix of H whose eigenvalues are identical to H. Thus, by the Cauchy interlacing theorem [START_REF] Elsner | A note on optimal block-scaling of matrices[END_REF],

D -1/2 HD -1/2 2 ≥ λ r+1 (H) = 1 + σ r+1 .
Applying similar arguments it can also be shown that

D 1/2 H -1 D 1/2 2 ≥ λ r+1 (H -1 ) = 1 1 -σ r+1 ,
which leads to the following lower bound on the spectral condition number,

cond 2 (R -T AR -1 ) ≥ 1 + σ r+1 1 -σ r+1 .
Finally, let us consider the specific choice

U 1 SV T 2 = R T 1 U r ΣV T r R 2 where U r and
V r are composed of the first r columns of U and V, respectively, in the SVD of

R -T 1 M R -1 2
. Consequently, we have σ i = σ i , i = 1, . . . , r. Thus, setting accordingly W = W and R = R we obtain by direct calculations:

cond 2 (R -T AR -1 ) = R -T AR -1 2 RA -1 R T 2 = 1 + σ r+1 1 -σ r+1 ,
and the proof is complete.

This manuscript is for review purposes only.

3. The Multilevel Weighted HODLR Preconditioning Scheme. In this section we introduce the multilevel HODLR preconditioning scheme for SPD matrices.

The method is based on the theory presented in section 2 and relaxation of the original problem. The motivation is twofold. First we demonstrate that the 1-level analysis can be extended to a multilevel preconditioning scheme. Second, we provide spectral bounds on the eigenvalues of the preconditioned system which give account for error propagation through the levels of the hierarchy.

In subsection 3.1 we give a brief review of the HODLR matrix structure which will be employed throughout the remainder of this paper. We focus on the symmetric case, since this work is concerned with the preconditioning of SPD matrices. In subsection 3.2 we introduce the preconditioning scheme, which is based on hierarchical construction and fast application of the inverse square roots, R -1 and R -T . In subsection 3.3 we briefly consider the associated memory and the computational costs of constructing and applying the scheme. An in-depth spectral analysis is presented in subsection 3.4. Our analysis provides estimates of the spectral bounds of the preconditioned system at each level, that take into account the approximation error at the lower levels of the hierarchy. These bounds are, however, of qualitative nature as they reflect a possible worst case scenario which is not likely to occur in practice.

A rigorous and detailed proof of the theory is given in subsection 3.5.

Symmetric HODLR Matrix Structure.

A symmetric HODLR matrix, K ∈ R n×n , can be described in the following recursive manner,

(17) K = K (0) 1 , K ( ) k = K ( +1) 2k-1 U ( +1) 2k-1 S ( ) k V ( +1) T 2k V ( +1) 2k S ( ) k U ( +1) T 2k-1 K ( +1) 2k ∈ R n ( ) k ×n ( )
k , for = 0, 1, . . . , L -1 and k = 1, 2, . . . , 2 , where is the level of K ( ) k in the hierarchy.

The triple products, U

( +1) 2k-1 S ( ) k V ( +1) T 2k
, represent low-rank blocks in the sense that ( 18)

U ( +1) 2k-1 ∈ R n ( +1) 2k-1 ×r ( ) k , V ( +1) 2k ∈ R n ( +1) 2k ×r ( ) k , S ( ) k ∈ R r ( ) k ×r ( ) k , where r
( ) k
is the rank of the corresponding off-diagonal block of K. Typically,

r ( ) k n ( +1) 2k-1 , n ( +1) 2k
. An illustration of the hierarchical structure of K is displayed in This manuscript is for review purposes only.

The common practice is to set the HODLR matrix as an approximation of a given matrix, A ∈ R n×n . The low-rank off-diagonal blocks of K satisfy For obtaining low-rank approximations satisfying [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF], the low-rank singular value decomposition (SVD) [START_REF] Golub | Matrix computations, Forth Edition[END_REF] which originated in [START_REF] Eckart | The approximation of one matrix by another of lower rank[END_REF] is, generally, considered the best choice, since it is optimal with respect to any unitarily invariant norm (2-norm, Frobenius). The computational cost required to obtain the SVD of M others, the rank revealing LU [START_REF] Miranian | Strong rank revealing LU factorizations[END_REF], rank revealing QR [START_REF] Gu | Efficient Algorithms for Computing a Strong Rank-Revealing QR Factorization[END_REF], randomized algorithms [START_REF] Frieze | Fast Monte-carlo Algorithms for Finding Low-rank Approximations[END_REF][START_REF] Liberty | Randomized algorithms for the low-rank approximation of matrices[END_REF][START_REF] Woolfe | A fast randomized algorithm for the approximation of matrices[END_REF], adaptive cross approximation [START_REF] Rjasanow | Adaptive cross approximation of dense matrices[END_REF] and boundary distance low-rank [START_REF] Aminfar | A fast block low-rank dense solver with applications to finite-element matrices[END_REF].

(19) M ( ) k -U ( +1) 2k-1 S ( ) k V ( +1) T 2k 2 ≤ τ ( ) k • M ( ) k
For more details see a review on this topic in [START_REF] Aminfar | A fast block low-rank dense solver with applications to finite-element matrices[END_REF].

3.2. Recursive Formula of the Multilevel Preconditioning Scheme. By our definitions [START_REF] Hackbusch | A Sparse Matrix Arithmetic Based on H-Matrices. Part I: Introduction to H-Matrices[END_REF] the principal blocks of K and the corresponding blocks of A are described by

(20) K ( ) k = K ( +1) 2k-1 U ( +1) 2k-1 S ( ) k V ( +1) T 2k V ( +1) 2k S ( ) k U ( +1) T 2k-1 K ( +1) 2k , A ( ) k = A ( +1) 2k-1 M ( ) k M ( ) T k A ( +1) 2k
, where K

( ) k , A ( ) k ∈ R n ( ) k ×n ( )
k , and the rank of each off-diagonal triple product approximation satisfies

(21) rk U ( +1) 2k-1 S ( ) k V ( +1) T 2k ≤ r ( ) k .
The key idea we propose is to construct each K ( )

k as an optimal 1-level preconditioning scheme of the matrix

B ( ) k = K ( +1) 2k-1 M ( ) k M ( ) T k K ( +1) 2k ∈ R n ( ) k ×n ( ) k ,
which is obtained by replacing the principal blocks of A ( )

k with the corresponding principal blocks of K ( )
k . This is a relaxation of the original problem, which facilitates a fast construction method. The resulting preconditioned global system condition number is no longer ensured to be minimal. However, the numerical results in section 4 indicate that the proposed approach is highly robust and, in general, attains superior condition number compared to HODLR approximations using other low-rank approximation schemes.

Before proceeding we introduce some necessary notations. Let R ( )

k ∈ R n ( ) k ×n ( ) k denote the square root of K ( ) k in the sense that K ( ) k = R ( ) T k R ( ) k . Let U ( ) k,r ∈
This manuscript is for review purposes only.

R n ( +1) 2k-1 ×r

( ) k and V ( ) k,r ∈ R n ( +1) 2k ×r ( ) k
be two thin matrices with orthogonal columns composed of the first r ( ) k left and right, respectively, singular vectors of the SVD, ( 22)

R ( +1) -T 2k-1 M ( ) k R ( +1) -1 2k = U ( ) k Σ ( ) k V ( ) T k , where Σ ( ) k,r ∈ R r ( ) k ×r ( ) k is the principal submatrix of Σ ( ) k = diag(σ ( ) k,1 , . . . , σ ( ) k,min{n ( +1) 2k-1 ,n ( +1) 2k }
) .

For brevity and clarity we will abuse the notation and employ U = U ( ) k,r and

V = V ( )
k,r . The proof of Theorem 1 implies that the following recursive formulas for a fast application of the inverse square roots hold,

R ( ) -T k = I + 1 2 U r U r V r -V r S + r -I 0 0 S - r -I U T r V T r U T r -V T r R ( +1) 2k-1 0 0 R ( +1) 2k -T , and
R ( ) -1 k = R ( +1) 2k-1 0 0 R ( +1) 2k -1 I + 1 2 U r U r V r -V r S + r -I 0 0 S - r -I U T r V T r U T r -V T r ,
where

S ± r = (I ± Σ ( ) k,r ) -1 2
. These formulas can be verified by writing the product

R ( ) -T k K ( ) k R ( ) -1 k
which, indeed, equals to the identity matrix, assuming K ( ) k is SPD. can be applied relatively fast in matrix product operations. In the case that a constant average rank, r ( ) = O(r), is taken for the off-diagonal blocks, the recursive representations implies the following relation ( 23)

C ( ) (m, n, r) = 2C ( +1) (m, n, r) + O m • n 2 • r ,
where C ( ) (m, n, r) denotes the computational cost of the operation R

( ) -T k • F ( ) k at level , F ( ) k ∈ R n ( )
k ×m , and n

( ) k = n ( ) = n/2 is assumed.
The first contribution in [START_REF] Parlett | The Symmetric Eigenvalue Problem[END_REF] stems from the recursive calls of the inverse square roots of the diagonal blocks.

The second contribution is associated with the cost of the matrix product operations.

Expanding [START_REF] Parlett | The Symmetric Eigenvalue Problem[END_REF] into a sum yields the total computational cost estimate,

C (0) (m, n, r) = log( n r ) =0 O (m • rn) = O (m • rn log (n)) ,
where the depth of the hierarchy is set by n = 2 L r. Similarly, the cost of storing

R (0) -1 k and R (0) -T k
is equal to O(rn log n) in the case where the average rank per level, r ( ) , is of O(r). See [START_REF] Aminfar | A fast block low-rank dense solver with applications to finite-element matrices[END_REF] for further details. As noted in the introduction if the HODLR scheme is generalized to HSS the costs to apply the preconditioner reduce to about O(n), see [START_REF] Xia | Effective and Robust Preconditioning of General SPD Matrices via Structured Incomplete Factorization[END_REF] for further details.

This manuscript is for review purposes only.

Constructing the preconditioner is accomplished by performing a single pass over the hierarchy from bottom to top. At each level we compute the low-rank factorizations of the triple products ( 24)

R ( +1) -T 2k-1 M ( ) k R ( +1) -1 2k
where k = 1, 2, . . . , 2 . Obtaining the low-rank factorization is performed by the following procedure:

• Capture the range of ( 24) in a matrix Q (k, ) L whose columns are orthonormal,

I -Q (k, ) L Q (k, ) T L • R ( +1) -T 2k-1 M ( ) k R ( +1) -1 2k ≈ 0 , Q (k, ) T L Q (k, ) L = I .
• Capture the range of the transpose of ( 24) in a matrix Q

(k, ) R
whose columns are orthonormal,

I -Q (k, ) R Q (k, ) T R • R ( +1) -T 2k M ( ) T k R ( +1) -1 2k-1 ≈ 0 , Q (k, ) T R Q (k, ) R = I .
• Compute the rank r (k, ) truncated SVD of the reduced matrix,

Q (k, ) T L • R ( +1) -T 2k-1 M ( ) k R ( +1) -1 2k • Q (k, ) R ≈ U ( ) k,r Σ ( ) k,r V ( ) k,r
• Reconstruct the WSVD left and right singular vectors matrices,

U ( ) k,r = Q L • U k,r , V ( ) k,r = Q R • V k,r .
If the effective rank of ( 24) is small, e.g. O(1), we can capture the range matrices,

Q (k, ) L and Q (k, ) R
, quickly by applying [START_REF] Rjasanow | Adaptive cross approximation of dense matrices[END_REF] and its traspose on a small set of randomized column vectors. See [START_REF] Martinsson | A Fast Randomized Algorithm for Computing a Hierarchically Semiseparable Representation of a Matrix[END_REF] for more details. However, if the effective rank of ( 24)

is not small, the procedure can become costly.

Spectral Estimates and Error Propagation.

Let us now focus on the spectral properties of the preconditioned submatrices, R

( ) -T k A ( ) k R ( ) -1 k
, where R

( ) k
is the square root of the principal submatrix K ( )

k . The submatrix A ( )
k is given in [START_REF] Liberty | Randomized algorithms for the low-rank approximation of matrices[END_REF]. Clearly the important case is = 0, since we are ultimately interested in preconditioning the input matrix, A = A

For brevity and clarity we will abuse the notation and employ the following representations in the spirit of section 2,

K ( ) k = K 1 U 1 SV T 2 V 2 SU T 1 K 2 , A ( ) k = A 1 M M T A 2 , R T i R i = K i (i = 1, 2) .
Note that R i represents an approximate square root of A i , as opposed to the exact square root that was assumed in section 2. We make the fundamental assumption that we have at our disposal spectral bounds estimates, (25)

α i ≤ λ min R -T i A i R -1 i ≤ λ max R -T i A i R -1 i ≤ β i (i = 1, 2) .
The lower-level bounds, α i and β i (i = 1, 2), can be obtained numerically, or possibly estimated analytically by the theory presented in this subsection. Note that in the case = L -1 we have α i = 1 = β i , and in the case = L -2 we have from section 2 the exact bounds

α i = λ min R -T i A i R -1 i = 1 -σ (L-1) 2k-2+i,r ∈ (0, 1] ,
This manuscript is for review purposes only.

β i = λ max R -T i A i R -1 i = 1 + σ (L-1) 2k-2+i,r ∈ [1, 2) .
The main result of the current subsection is presented in Theorem 2. The theorem provides a description of the behavior of the current-level spectral bounds, ( 26)

α ≤ λ min R ( ) -T k A ( ) k R ( ) -1 k ≤ λ max R ( ) -T k A ( ) k R ( ) -1 k ≤ β ,
as a function of the lower-level bounds [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF] and the rank of the off-diagonal blocks,

r = r ( )
k . The definition of the bounds α and β (26) is based on variational formulation and provided in Lemma 2. The analysis requires sufficient (but not necessary) conditions on the given lower-level bounds, α i and β i (i = 1, 2). We show that the proposed HODLR preconditioning scheme, essentially, maps both the r largest and the r smallest eigenvalues to a closed segment containing 1. When this segment is small, the preconditioner retains optimality or near optimality. We also show that the sensitivity of the spectral bounds to the inaccuracies

K i = A i (i = 1, 2) is governed
by the Cauchy-Bunyakowski-Schwarz (CBS) constant [START_REF] Axelsson | Algebraic multilevel preconditioning methods. I[END_REF][START_REF] Axelsson | Algebraic Multilevel Preconditioning Methods, II[END_REF].

Theorem 2. Let A = A 1 M M T A 2 , K = K 1 U 1 SV T 2 V 2 SU T 1 K 2 ,
be symmetric matrices of the same dimensions and partition where A is SPD. Assume that the off-diagonal triple product approximation U 1 SV T 2 satisfy

U 1 = R T 1 U r , S = Σ r , V 2 = R T 2 V r ,
where U r and V r are composed of the first r left and right, respectively, singular vectors of the SVD,

R -T 1 M R -1 2 = UΣV T , Σ = diag(σ 1 , . . . , σ min{n1,n2} ) , and
R T i R i = K i (i = 1, 2).
Assuming there exist real positive constants,

(27) 0 < α 1 , α 2 ≤ 1 ≤ β 1 , β 2 ,
such that

(28) 0 < α i x T i K i x i ≤ x T i A i x i ≤ β i x T i K i x i ∀ x i ∈ R ni ,
we have the following spectral estimates:

1. If σ 1 < √ α 1 α 2 then (29) α = min α avg 1,2 -σ 2 r+1 + α dif 1,2 2 , α avg 1,2 + σ 1 1 + σ 1 -δ 1,2 α , where
α avg 1,2 = α 1 + α 2 2 , α dif 1,2 = α 1 -α 2 2 , (α = α or β) , and
δ 1,2 α = 1 4 α avg 1,2 -σ 1 1 -σ 1 - α avg 1,2 + σ 1 1 + σ 1 2 + α dif 1,2 2 1 -σ 2 1 - 1 2 α avg 1,2 -σ 1 1 -σ 1 - α avg 1,2 + σ 1 1 + σ 1 ,
is a positive lower spectral bound of the preconditioned system,

0 < α ≤ λ R -T AR -1 , K = R T R .
This manuscript is for review purposes only.

If σ

1 < 1 then (30) β = max β avg 1,2 + σ 2 r+1 + β dif 1,2 2 , β avg 1,2 -σ 1 1 -σ 1 + δ 1,2 β ,
where

β avg 1,2 = β 1 + β 2 2 , β dif 1,2 = β 1 -β 2 2 , (β = α or β) ,
and

δ 1,2 β = 1 4 β avg 1,2 -σ 1 1 -σ 1 - β avg 1,2 + σ 1 1 + σ 1 2 + β dif 1,2 2 1 -σ 2 1 - 1 2
β avg 1,2 -σ 1 1 -σ 1 - β avg 1,2 + σ 1 1 + σ 1 .
is an upper spectral bound of the preconditioned system,

0 < λ R -T AR -1 ≤ β , K = R T R .
Remarks. The justification for ( 27) is a consequence of Theorem 2, which shows that α (29) and β (30) are monotonically non-increasing and non-decreasing as functions of the level, respectively. This observation is also supported by numerical evidence in section 4. If √ α 1 α 2 ≤ σ 1 < 1, then the preconditioned system remains SPD.

However, the theory presented here can not predict the positive value of the lower spectral bound, α (29).

From Theorem 2 we observe that each estimated bound, α or β, is a minimum or a maximum, respectively, of two competing terms: the first depends on the largest singular value, σ 1 , and the second is a function of the truncation error, σ r+1 . In fact, when the truncation error becomes sufficiently small it does not affect the values of the bounds, which are governed solely by the terms depending on the largest singular value. Thus in this case, improving the approximation by increasing the rank r does not improve the corresponding condition number estimate, β/α. An illustration of this observation is given in Figure 3.

The last observation as displayed in Figure 3 indicates that the value of σ 1 is central to the estimation of the spectral bounds, and effectively dominates the condition number of the preconditioned system. In this sense σ 1 reflects the sensitivity of the condition number of R

( ) -T k A ( ) k R ( ) -1 k
to the lower level inaccuracies. It can be shown that σ 1 is the so-called Cauchy-Bunyakowski-Schwarz (CBS) constants of the matrix K, which is defined by ( 31)

σ 1 = sup x1,x2 =0 x T 1 M x 2 x T 1 K 1 x 1 x T 2 K 2 x 2 ≥ 0 .
Definition (31) originated from the theory of Algebraic Multilevel Iterations Methods [START_REF] Axelsson | Algebraic multilevel preconditioning methods. I[END_REF][START_REF] Axelsson | Algebraic Multilevel Preconditioning Methods, II[END_REF], and coincides with the principal angle (cosine of the smallest angle) between the column space of I1 0 T and the column space of 0 I2 T with respect to the inner product x, y A = y T Ax. Thus, σ 1 represents the local contribution of the upper level to the overall condition number. Using [START_REF] Yang | Sparse Hierarchical Solvers with Guaranteed Convergence[END_REF] with assumption (27) leads to the following relation

1 √ β 1 β 2 ≤ σ 1 σ exact 1 ≤ 1 √ α 1 α 2 ,
where σ exact 1.0

β avg 1,2 -σ 1 1-σ 1 1.5 1.75
Relative Truncation Error:

σ r+1 /σ 1 Spectral Bounds α avg 1,2 -σ1 β avg 1,2 + σ1
Fig. 3: Spectral Bounds. A typical behavior of the spectral bounds displayed for the case α 1 = α 2 and β 1 = β 2 . The lower bound α (29) vs. σ r+1 /σ 1 is plotted in blue, and the upper bound β (30) vs. σ r+1 /σ 1 is plotted in red.

given matrix, A, and the chosen partition. If K is close to A then we can expect σ 1 to be close to σ exact 1 , and in this case we have little influence over its value.

Regarding the spectrum of the preconditioned system, the interpretation of Theorem 2 is similar to the interpretation of Theorem 1. From the proof it can be inferred that two-sided block Jacobi (i.e., the case r = 0) effectively maps the spectra of the bounding preconditioned systems to two segments centered around α avg 1,2 and β avg 1,2 ,

α avg 1,2 -σ 2 1 + α dif 1,2 2 , α avg 1,2 + σ 2 1 + α dif 1,2 2 , β avg 1,2 -σ 2 1 + β dif 1,2 2 , β avg 1,2 + σ 2 1 + β dif 1,2 2 .
The multilevel Weighted HODLR preconditioning scheme does the same, but also maps the largest and smallest eigenvalues to the segments (32)

σ 1 + α avg 1,2 1 + σ 1 -δ 1,2 α , σ 1 -α avg 1,2 1 -σ 1 + δ 1,2 α , (33)
β avg 1,2 + σ 1 1 + σ 1 -δ 1,2 β , β avg 1,2 -σ 1 1 -σ 1 + δ 1,2 β ,
respectively. Thus, assuming the segments (32) and (33) are small, a significant improvement in the condition number as well as the clustering of the spectrum of the original preconditioned system is expected. An illustration is given in Figure 4. The figure is similar to Figure 1 where the main difference is that the weighted HODLR preconditioning scheme now maps the extreme eigenvalues to an interval containing 1 and not exactly to 1.

This manuscript is for review purposes only.

A BJ Weighted

λ min α BJ α avg 1,2 -σ r+1
1.0

β avg 1,2 + σ r+1

Approximate Bounds

Fig. 4: Spectrum Clustering for the Multilevel Weighted HODLR Preconditioning Scheme. The spectrum of some SPD matrix A and the transformation it goes after preconditioning by block Jacobi (BJ) and the multilevel weighted HODLR preconditioning scheme are displayed. The spectra are ordered from the left to the right starting from A, followed by Block Jacobi (BJ) and ends up with the multilevel weighted HODLR scheme. The spectral bounds α (29) and β (30) are marked on the right y-axis, while the spectral bounds for the block Jacobi case α BJ = α(r = 0) and β BJ = β(r = 0) are marked on the left y-axis.

Lemma 2. Let

A = A 1 M M T A 2 , K = K 1 U 1 SV T 2 V 2 SU T 1 K 2 ,
and the associated lower-level bounds

(34) 0 < α 1 , α 2 ≤ 1 ≤ β 1 , β 2 ,
satisfy the assumptions of Theorem 2.

Let us define

(35) K = α 1 K 1 M M T α 2 K 2 , K = β 1 K 1 M M T β 2 K 2 .
Then we have:

1. The matrices K, K are SPD iff

σ 1 < √ α 1 α 2 , σ 1 < β 1 β 2 ,
respectively, where σ 1 is the largest singular of R -T 1 M R -1 2 .

2. If K, K are SPD, there exist two positive constants, α and β, such that

(36) α = min x =0 x T Kx x T Kx ≤ λ min ( B T A B) ≤ λ max ( B T A B) ≤ max x =0 x T Kx x T Kx = β .
Proof. of Lemma 2.

To show the first part of the lemma we consider a general matrix of the form

H = δ 1 K 1 M M T δ 2 K 2 , δ 1 , δ 2 > 0 .
Let us apply the following two-sided transformations

H = 1 √ δ1 R -T 1 0 0 1 √ δ2 R -T 2 H 1 √ δ1 R -1 1 0 0 1 √ δ2 R -1 2 = I 1 √ δ1δ2 M 1 √ δ1δ2 M T I ,
where

M = R -T 1 M R -1 2 .
The matrix H is SPD iff H is SPD as well. Thus, by Lemma 1 the matrix H is SPD iff 1 -σ 1 / √ δ 1 δ 2 > 0, and the conditions ensuring K, K, and K are SPD immediately follow.

For the second part of the lemma it is sufficient to assume that K is SPD which, by the first part, ensures that K and K are SPD as well. Accordingly, we obtain the following inequalities

x T Kx x T Kx ≤ x T Kx x T Ax ≤ x T Kx x T Kx ∀ x = 0 .
The Lagrangian stationary points of each generalized Rayleigh quotient in the inequalities above constitute the spectrum of each preconditioned system. Thus, the proof is complete.

Lemma 3. Let H = D (1) D (2) D (2) D (3) ∈ R 2r×2r
, where D (i) (i = 1, 2, 3) are diagonal matrices,

D (i) = diag(d (i) 1 , . . . , d (i) r ) .

If d

(2) j = 0 for all j = 1, 2, . . . , r, then

spec(H) = λ - j m j=1 ∪ λ + j m j=1 , λ ± j = d (1) j + d (3) j 2 ± d (1) j -d (3) j 2 2 + d (2) j 2 ,
where spec(H) denotes the spectrum of the symmetric matrix H.

Proof. of Lemma 3.

From the given structure of H it is clear that λ ∈ R is an eigenvalue of H iff for some j = 1, 2, . . . , m the vectors (d

(1) j -λ, d (2)
j ) and (d

(2) j , d (3)
j -λ) are linearly dependent.

Since we have assumed d

(2) j = 0, we have that (d

(1) j -λ, d (2)
j ) and (d

(2) j , d (3)
j -λ) are linearly dependent iff d (1) j -λ d (2) j = d (2) j d (3) j -λ ⇔ (d (1)
j -λ)(d (3)
j -λ) -(d (2) j ) 2 = 0 .
The solution to the quadratic equation above is

λ = λ ± j = d (1) j + d (3) j 2 ± d (1) j -d (3) j 2 2 + d (2) j 2 ,
and the proof is complete.

This manuscript is for review purposes only.

Proof. of Theorem 2.

Considering the conditions of Theorem 2 we have by Lemma 2 that the spectral bounds, α and β, satisfying 0 < α = min

x =0

x T Kx x T Kx ≤ λ min (R -T AR -1 ) ≤ λ max (R -T AR -1 ) ≤ max

x =0

x T Kx x T Kx = β , exist where K and K are defined in Lemma 2.

To find the exact values of α and β, we consider a generalized Rayleigh quotient

Q(x) = x T Hx x T Kx , H = δ 1 K 1 M M T δ 2 K 2 ,
whose range is strictly positive. Thus, Q(x) represents either x T Kx/x T Kx or x T Kx.

We apply the change of variables, x =

R -T 1 0 0 R -1 2
ξ, and obtain the following equivalent representation

Q(x) = ξ T Hξ ξ T Kξ , H = δ 1 I 1 M M T δ 2 I 2 , K = I 1 M r M T r I 2 ,
where

M = R -T 1 M R -1 2 and M r = U r Σ r V T r is the r-rank weighted SVD approxima- tion of M.
Let w i denote the i-th column of the orthogonal matrix W (13) as defined in Lemma 1. Then we have:

1. Kw i = (1 + σ 1 )w i , i = 1, 2, . . . , r.
2. Kw n1+i = (1 -σ 1 )w n1+i , i = 1, 2, . . . , r.

3. Kw j = w j , j = 1, . . . , r, n 1 + 1, . . . , n 1 + r.

and similarly for H, it can be verified that:

1. Hw i = (δ avg 1,2 + σ 1 )w i + δ dif 1,2 w n1+i , i = 1, 2, . . . , min{n 1 , n 2 }.

2. Hw n1+i = (δ avg 1,2 -σ 1 )w n1+i + δ dif 1,2 w i , i = 1, 2, . . . , min{n 1 , n 2 }.

3. Hw j = w j , j = 1, . . . , min{n 1 , n 2 }, n 1 + 1, . . . , n 1 + min{n 1 , n 2 }.

where δ avg 1,2 = (δ 1 + δ 2 )/2 and δ dif 1,2 = (δ 1 -δ 2 )/2. Clearly, both K and H are invariant over the subspaces Z = span{w 1 , . . . , w r , w n1+1 , . . . , w n1+r } and its orthogonal complement, Z ⊥ . Hence, by the properties of the generalized Rayleigh quotient we have:

max x =0 Q(x) = max max ξ∈Z\{0} Q(x) , max ξ∈Z ⊥ \{0} Q(x) ,
and

min x =0 Q(x) = min min ξ∈Z\{0} Q(x) , min ξ∈Z ⊥ \{0} Q(x) .
By our results so far, if x = ξ ∈ Z ⊥ then Q(x) = ξ T Hξ/ξ T ξ. Let us apply the change of variables of the form ξ = Cζ ∈ Z ⊥ , given explicitly by

ξ = ζ 1 w r+1 + . . . + ζ n1-r w n1 + ζ n1-r+1 w n1+r+1 + . . . + ζ n1+n2-2r w n1+n2 ,
where ζ i is the i-th coordinate of ζ and as before w i denotes the i-th column in the orthogonal matrix W (13). Then, for any ξ ∈ Z ⊥ we obtain

Q(x) = ζ T H Z ⊥ ζ ζ T ζ , H Z ⊥ = D (1) Z ⊥ D (2) Z ⊥ D (2) T Z ⊥ D (3) Z ⊥ ,
This manuscript is for review purposes only.

where D

(2)

Z ⊥ = δ dif 1,2 I n1,n2 and
D ( 1
)
Z ⊥ = diag δ avg 1,2 + σ r+1 , . . . , δ avg 1,2 + σ n1 if n 1 ≤ n 2 , diag δ avg 1,2 + σ r+1 , . . . , δ avg 1,2 + σ n1 , δ 1 , . . . , δ 1 if n 1 > n 2 . D (3)
Z ⊥ = diag δ avg 1,2 -σ r+1 , . . . , δ avg 1,2 -σ n2 if n 2 ≤ n 1 , diag δ avg 1,2 -σ r+1 , . . . , δ avg 1,2 -σ n2 , δ 2 , . . . , δ 2 if n 2 > n 1 .
Now,by Lemma 3, we obtain that the spectrum of H Z ⊥ contains the sets

δ avg 1,2 + σ 2 r+1 + δ dif 1,2 2 , . . . , δ avg 1,2 + σ 2 min{n1,n2} + δ dif 1,2 2 , , δ avg 1,2 -σ 2 min{n1,n2} + δ dif 1,2 2 , . . . , δ avg 1,2 -σ 2 r+1 + δ dif 1,2 2
.

Hence, we conclude that min

ξ∈Z ⊥ \{0} Q(x) = δ avg 1,2 -σ 2 r+1 + δ dif 1,2 2 , max ξ∈Z ⊥ \{0} Q(x) = δ avg 1,2 + σ 2 r+1 + δ dif 1,2 2 .
For the case ξ ∈ Z let us apply the change of variables of the form ξ = Cψ ∈ Z, given explicitly by ξ = ψ 1 w 1 + . . . + ψ r w r + ψ r+1 w n1+1 + . . . + ψ 2r w n1+r , where ψ i is the i-th coordinate of ψ and as before w i denotes the i-th column in the orthogonal matrix W (13). Then, for any ξ ∈ Z we obtain

Q(x) = ψ T H Z ψ ψ T ψ , H Z = D (1) Z D (2) Z D (2) Z D (3) Z , D (1)
Z = diag δ avg 1,2 + σ 1 1 + σ 1 , . . . , δ avg 1,2 + σ r 1 + σ r , D (2)
Z = diag δ dif 1,2 1 -σ 2 1 , . . . , δ dif 1,2 1 -σ 2 r , D (3)
Z = diag δ avg 1,2 -σ 1 1 -σ 1 , . . . , δ avg 1,2 -σ r 1 -σ r .
Applying once more the outcome of Lemma 3 we have that the spectrum of H Z is composed of the following values 1 2

δ avg 1,2 + σ i 1 + σ i + δ avg 1,2 -σ i 1 -σ i ± 1 4 δ avg 1,2 + σ i 1 + σ i + δ avg 1,2 -σ i 1 -σ i 2 + δ dif 1,2 2 1 -σ 2 i
, where i = 1, 2, . . . , r and the proof is complete.

This manuscript is for review purposes only.

4. Numerical Study. This section contains the experimental part of this work.

The main goal is to demonstrate the effect of different low-rank approximations [START_REF] Hackbusch | A Sparse H-Matrix Arithmetic. Part II: Application to Multi-Dimensional Problems[END_REF] for the off-diagonal blocks on the preconditioned system using HODLR. We perform a comparative study and consider the following low-rank techniques:

• R-HODLR: the off-diagonal low-rank factorizations are obtained in the standard or regular approach using truncated SVD.

• C-HODLR: the off-diagonal low-rank factorizations are obtained using truncated SVD with additional imposed constraints as described in [START_REF] Bebendorf | Hierarchical matrix approximation with blockwise constraints[END_REF].

• W-HODLR: the off-diagonal low-rank factorizations are obtained using the weighted HODLR preconditioning scheme for the multi-level case. The construction and application of the scheme follows the outlined procedure in subsection 3.3.

Employing SVD is done for convenience and uniformity of the comparison, and can be replaced, in practice, by other more efficient low-rank approximation techniques.

Subsection 4.1 describes the computational setting, and presents a pair of severely ill-conditioned sparse systems which have been used in the numerical simulations. Another simplified numerical example along a more detailed analysis can be found in [START_REF] Agullo | Nearly optimal fast preconditioning of symmetric positive definite matrices[END_REF].

In subsection 4.2 we describe the numerical results of the PCG approximation using HODLR preconditioning. The results indicate that the weighted low-rank factorization scheme proves to be superior to other standard techniques.

4.1. Sparse Matrices and Computational Setting. In the presented numerical study we have explored and analyzed the PCG solution for the following sparse matrices, which have been picked from the SuiteSparse matrix collection [START_REF] Davis | The University of Florida Sparse Matrix Collection[END_REF].

• bcsstk16: 4, 884 × 4, 884 , SPD , spectral condition number ≈ 4.94 • 10 9 .

• bcsstk15: 3, 948 × 3, 948 , SPD , spectral condition number ≈ 6.53 • 10 9 .

For constructing the HODLR preconditioning schemes we interpret each matrix as a discrete graph and apply a balanced partitioning using successive bisections. We employ Scotch [START_REF] Chevalier | PT-Scotch: A tool for efficient parallel graph ordering[END_REF] for each bisection dividing a given vertex set into two distinct sets of approximately equal size whose cut is minimal, i.e., the number of edges running between the separated subsets is as small as possible. The process starts with the entire set of vertices, and then applied recursively on each separated subset until reaching the predetermined bottom level of the hierarchy, L.

Construction of the preconditioning schemes and the iterative solution of the preconditioned system has been implemented with a Fortran90 code. In all the simulations we have employed the following selections:

• L = log 2 (n/100) for an n × n matrix as the lowest level of the hierarchy, which forces the size of the smallest blocks in the partition under 100.

• Constant off-diagonal block ranks over all levels of the hierarchy with the following O(1) values:

(37) r [START_REF] Liberty | Randomized algorithms for the low-rank approximation of matrices[END_REF], and then computing the low-rank factorization [START_REF] Rjasanow | Adaptive cross approximation of dense matrices[END_REF] on the reduced block. For the sparse case, this procedure is, typically, equivalent in terms of complexity to the randomized technique [START_REF] Martinsson | A Fast Randomized Algorithm for Computing a Hierarchically Semiseparable Representation of a Matrix[END_REF].

This manuscript is for review purposes only. where i = 1, 2, . . . is the iteration step index and x (i) denotes the approximate solution at step i. The results indicate that in all the test cases, the W-HODLR scheme outperforms the other techniques, and retain good properties even when low accuracy for the off-diagonal blocks approximations is employed.

Figure 5 contains plots of the PCG convergence history profiles for bcsstk16.

All plots in this case indicate that increasing the constant rank (37), improves the approximation quality, and achieves faster convergence rate. It is also evident that R-HODLR and C-HODLR achieve similar convergence with the same memory resources, while W-HODLR converges faster with the same memory resources. Figure 6 contains plots of the PCG convergence history profiles for bcsstk15. The results show that R-HODLR and C-HODLR fail to converge in 1000 PCG iterations. In fact, setting constant rank 0, i.e., Block Jacobi, performs better then using these schemes with a constant rank greater than zero. This occurs because the use of naive approximations for the off-diagonal blocks makes the problem even more ill-conditioned. The W-HODLR scheme, however, converges with excellent convergence rates, where the convergence rate improves when the constant rank (37) is increased. We have presented a mathematical analysis for obtaining low-rank factorizations, that minimize the spectral condition number of the preconditioned system for the 1level (2 × 2) case. The key idea was to properly reweight the blocks prior the low-rank factorization, which leads to a minimum spectral condition number.

The presented theory has been extended to HODLR preconditioning schemes, including analysis of the spectral properties and bounds that take into account the error propagation through the levels of the hierarchy.

The numerical experiments indicate, that employing proper reweighting for the off-diagonal blocks prior to low-rank compression, leads to a HODLR preconditioning scheme that requires far less computational resources for the same quality of performance of convergence than using the other low-rank compression techniques.

As noted in the introduction a major goal of this work is to provide an analysis of optimal choice of low-rank approximations for a simple case; i.e., HODLR, which could lead to an extended analysis for the strong hierarchical case. This point will be explored and pursued in a future study. This manuscript is for review purposes only.

Fig. 1 :

 1 Fig. 1: Spectrum Clustering of the Optimal 1-level Preconditioning Scheme.The spectrum of some SPD matrix A and the transformations it goes after preconditioning by block Jacobi (BJ) and the optimal 1-level preconditioning scheme are displayed. The spectra are ordered from the left to the right starting from A, followed by BJ and end up with the optimal scheme.

Figure 2 Fig. 2 :

 22 Figure 2.

2 ,k

 2 denotes the corresponding off-diagonal block of A and τ ( ) k > 0 is a chosen tolerance. Typically, a prior reordering of the matrix rows and columns is needed to confirm that r ( ) k are, indeed, low.

kk / 2 .

 2 is relatively expensive necessitating an O m 3 operations, where m = n ( ) For this reason a variety of fast approximation algorithms attempting to efficiently obtain a low-rank approximation close enough to the low-rank SVD have been proposed. These include, among

3. 3 .

 3 Utilization and Construction of the Preconditioning Scheme. The recursive representations of R

1 1 are correlated where σ exact 1 is

 11 is the corresponding CBS constant of A. The important conclusion here is that σ 1 and σ exact intrinsically predetermined by the This manuscript is for review purposes only.

3. 5 .

 5 Proof of Theorem 2. The proof of Theorem 2 is based on Lemma 2 which provides the definition of the bounds α (29) and β (30), and on the technical result presented in Lemma 3.

1 , 2 .

 12 . . . , L , k = 1, . . . , Note that r ( ) k ≡ 0 reduces the preconditioning scheme to block Jacobi (BJ), regardless of the specific low-rank factorization technique. The construction of the low-rank factorizations (24) follows the path outlined in subsection 3.3. We have produced fast low-rank factorizations by first removing all the zero rows and columns of the sparse block M ( ) k

4. 2 .

 2 Numerical Results and Analysis. This subsection contains the numerical results of PCG solution for the chosen sparse systems, using R-HODLR, C-HODLR and W-HODLR preconditioning schemes. For both matrices, bcsstk16 and bcsstk15, we have set the right-hand side to b = 1, and the iterative approximation was stopped at the first occurrence ofAx (i) -b 2 ≤ 10 -8 b 2 ,

5 .

 5 Summary and Future Work. In this work we have addressed the problem of choosing low-rank factorizations in fast hierarchical algorithms for preconditioning SPD matrices.

25 Fig. 5 :

 255 Fig.5: PCG Convergence History: 'bcsstk16'. Three plots showing PCG convergence history profiles, i.e., the values Ax (i) -b 2 / b 2 as a function of the iteration number i, for each preconditioning scheme. Each plot displays various profiles, where each profile corresponds to a different constant rank value (37) of the approximations for the off-diagonals blocks by low-rank factorizations.

25 Fig. 6 :

 256 Fig.6: PCG Convergence History: 'bcsstk15'. Three plots showing PCG convergence history profiles, i.e., the values Ax (i) -b 2 / b 2 as a function of the iteration number i, for each preconditioning scheme. Each plot displays various profiles, where each profile corresponds to a different constant rank value (37) of the approximations for the off-diagonals blocks by low-rank factorizations.

This manuscript is for review purposes only.

E. AGULLO, E. DARVE, L. GIRAUD, AND Y. HARNESS

Acknowledgments. We would like to thank the anonymous referees for their valuable remarks, questions and comments that enabled us to substantially improve this paper.

This manuscript is for review purposes only.