Yifat Cohen

Moran Levi

Uri Lesmes

Marielle Margier

Emmanuelle Reboul

Yoav D Livney

Re-assembled casein micelles improve in vitro bioavailability of vitamin D in a Caco-2 cell model

The pandemic of vitamin D (VD) deficiency, and the global rise in obesity stimulate a need for staple low-fat foods and beverages enriched with VD. In light of consumer demand for a clean label, the use of natural endogenous food ingredients as delivery vehicles is of great interest. To this end, reassembled casein micelles (rCM) have been shown to help retain VD during processing and shelf life as well as provide high bioavailability in low-fat milk and non-fat yoghurt. This follow up study focused on the performance of VD-loaded rCM after drying and reconstitution, considering VD retention during simulated digestion, and the subsequent in-vitro bioavailability of the vitamin. RCM conferred great protection to VD 3 during simulated digestion with significant (p<0.05) increase in vitamin retention for 1 h under gastric conditions. This observation is belived to arise from the vitamin-casein binding and the system's natural gelation (curd formation) near the casein isoelectric point that seclude the vitamin from environmental stressors and couple its release with digestive proteolysis of the rCM matrix. Vitamin absorption by Caco-2 cells from digested rCM was not significantly different from absorption of digested free VD.

However, thanks to the highly protective effect of the rCM, against VD gastric degradation, the overall effect of the rCM was a 4-fold higher bioavailability, compared to the free VD.

Introduction

Increasing evidence support the linkage between diet and human health 1 . This has led to numerous attempts to develop functional foods that carry added nutritional value. Such attempts include the enrichment and fortification of food with bioactive ingredients such as vitamins, polyphenols, essential fatty acids and carotenoids 2 . However, the incorporation of bioactive ingredients in food and drink systems poses several challenges. First, the enrichment of aqueous food or drink systems harbours a challenge of dosing accuracy and uniformity to avoid overdose, mainly in the case of oil-soluble bioactives, due to their poor solubility. Moreover, the incorporation of a bioactive ingredient within the food or drink system may result in its degradation and loss of activity during processing, product shelf life and digestion. Bioactive ingredients may cause unwanted effects on the sensory attributes of the enriched food system, and they can also react with other components present in the enriched system, which may decrease (or increase) their bioavailability [1][2][3] . To overcome these challenges, various bioactive delivery systems are being developed 2 . Delivery systems are commonly classified according to the main building blocks, which are necessarily GRAS (generally recognized as safe) food materials, preferably natural, and ideally-endogenous components of the product, thereby shortening the ingredient list and minimizing changes to the product sensory attributes, for higher consumer satisfaction. There are four main groups of delivery systems in food: protein-based, carbohydrate-based, lipid-based, and mixed systems 2 . Due to the negative health implications of high fat consumption 4,5 , there is a need for enriching low-fat and non-fat food and drinks with health-promoting oil-soluble bioactives, preferably using non lipid, amphiphilic delivery systems.

Proteins are natural and diverse class of macronutrients, with various origins and techno-functionalities.

Therefore, several studies have examined the development and use of protein-based delivery systems for lipophilic nutraceuticals [6][7][8][9][10][11][12][13] . Pertinent to vitamin D (VD), soybean β-conglycinin 10 , corn protein hydrolysate 11 , hydrophobins 9 , potato-protein 14 and milk proteins 8,12 have all shown good performances for the delivery of VD. Out of these, milk proteins seem particularly effective tools to deliver vitamins, minerals and lipophilic compounds, due to their evolutionary role in milk, which provides total nutrition to new-borns 15 .

The recommended dietary allowance (RDA) of VD is 600 international units (IU)/day for children over 12 months of age, adolescents and adults according to Institute of medicine report 16 . The RDA for elderly people above 70 years is 800 IU/day 16 . Nowadays, the actual daily consumption of VD is low and does not exceed 400 international units (IU), in most populations 17,[START_REF] Kiely | Scandinavian journal of clinical and laboratory investigation[END_REF] . Modern lifestyle and minimization of exposure to sun light render natural production of VD in the skin (by UVB radiation) to be generally insufficient, even in sunny countries 19 . VD is not only biologically important but also associated with prevention of various diseases 20,21 , therefore, its supplementation to foods is a sensible public health measure.

This study focuses on the potential use of caseins, and particularly of re-assembled casein micelles (rCM), as delivery vehicles for VD and extend a few recent studies. First, rCM were shown to be suitable for the enrichment of reduced-fat products with VD and other oil-soluble nutraceuticals 8,22,23 . RCM were found to protect VD 3 during UV exposure 22 , thermal treatment (80°C, 1 min) and during 28 days of cold storage 8 . Then, a clinical human study concluded that the bioavailability of VD encapsulated in rCM within 1% fat milk was at least as good as that of VD in a commercial dietary supplement aqueous formulation in which VD 3 was emulsified using Tween-80 8 . In another clinical trial, the enrichment of fatfree yoghurt with VD 3 entrapped in rCM showed similarly high bioavailability as when using the synthetic surfactant Tween-80 24 , though significantly better sensory properties were obtained with rCM. Another recent clinical study found that the bioavailability of VD in rCM in non-fat yoghurt was not significantly different from its bioavailability within the fat (in a 3% fat yoghurt), or in its presence (VD in rCM, in 3% fat yoghurt) 25 .

The objectives of this follow up study were to evaluate the performance of the VD-loaded rCM after drying and reconstitution, and to gain deeper understanding of their behaviour during simulated digestion and of the consequent in-vitro bioavailability of the VD encapsulated in rCM.

Materials and Methods

Materials

Commercial sodium caseinate (Casinella-QN 94.2% protein Molkerei Meggle Wasserburg GMBH & Co. KG) was kindly donated by Kelta Ltd. Israel and Molkerei Meggle Wasserburg GMBH & Co. KG, Germany. VD 3 in powder (USP grade), potassium phosphate dibasic (USP), potassium citrate tribasic monohydrate (FCC) and calcium chloride dihydrate (FCC), pepsin (Sigma P7000, ≥250 units/mg protein), Trypsin (Sigma T0303, 13,000-20,0000 BAEE units/mg protein), chymotrypsin (Sigma C4129, ≥40 units/mg protein), pancreatic lipase (Sigma L3126, 100-500 units/mg protein), sodium glycodeoxycholate (Sigma G9910) and taurocholic acid salt hydrate (Sigma T4009) were purchased from Sigma-Aldrich (Rehovot, Israel). Sodium chloride and Sodium bicarbonate were purchased from Frutarom Industries Ltd. (Herzliya, Israel) and Hydrochloric acid was purchased from Gadot group (Netanya, Israel). Dulbecco's modified Eagle's medium (DMEM) containing 4.5 g/L glucose, penicillin/streptomycin, trypsin-EDTA (500 mg/L and 200 mg/L, respectively) and PBS were purchased from Life Technologies (Saint Aubin, France). Fetal bovine serum was obtained from PAA (Vélizy Villacourblay, France).

Methods

VD 3 -casein binding characterization, VD 3 -loaded rCM preparation, rCM freeze-drying and reconstitution

Casein solution and VD 3 loading: A solution of 6.11 mg/ml sodium caseinate was prepared with double distilled water and stirred overnight at room temperature for complete hydration prior to vitamin loading. Samples were prepared by adding ethanolic VD 3 solution (8 mg/ml) into a vial containing either double distilled water or casein solution while vortexing, followed by gently stirring overnight at 4°C. Final VD 3 and caseinate concentrations in the samples were 260 µM each. A control of casein solution at the same concentration without VD 3 was prepared by the same procedure except that an equivalent amount of pure ethanol was added to the pre-dissolved sodium caseinate. Final ethanol concentration in the samples was 1.2%.

Light microscopy: Light microscopy images of the pure VD 3 dispersions in water (containing 1.2% ethanol) were compared to those of VD 3 : casein 1:1 molar ratio solutions. Images were taken at ×40 magnification using Olympus BX51 light microscope, operated in bright-field optical mode.

Binding studies by UV absorbance spectra measurements: The binding studies were performed according to the methods described by Israeli-Lev et al. 9 . UV-vis spectra were measured by UV spectrophotometer (Ultrospec 3000, GE Healthcare). All samples were diluted 10-fold before measurements. Measurements were made in duplicate at 23°C. Preparation of VD 3 loaded rCM: VD 3 loaded rCM were prepared according to the protocol described by Haham et al. 8 Briefly, to the solution of sodium caseinate loaded with VD 3 , described above, K 2 HPO 4 and tri potassium citrate solutions were added, followed by drop wise addition of CaCl 2 while stirring. Then, the samples were homogenized using a Micro-DeBee (Bee Int'l Inc., South Easton, MA, USA) ultra-high pressure homogenizer.

Freeze drying and reconstitution: 5% [wt/wt] Maltodextrin (MD, DE-18) was added to the rCM dispersion, as a cryoprotectant/drying aid, followed by stirring for 30 min at room temperature.

Suspensions were then quench-frozen in liquid nitrogen and freeze-dried. Reconstitution was achieved by adding distilled water and stirring overnight at 4°C. Particle size distribution of the rCM before and after freeze-drying and reconstitution were determined by dynamic-light-scattering (DLS) using NICOMP TM 380, Particle size/zeta potential analyser (PSS, Santa Barbara, CA, USA). Measurements were made at 23°C. Particle size distributions were analysed by the NICOMP TM volume-weighted distribution analysis. VD 3 concentration before and after freeze-drying and reconstitution was evaluated by reversed phase (RP) HPLC.

Evaluation of VD 3 bioaccessibility using an in vitro simulated digestion model

Simulated digestion process: To assess casein digestibility and VD 3 retention a semi-dynamic in vitro digestion model was applied and comprised of two steps: first a 1h semi-dynamic gastric digestion followed by a two hour intestinal digestion step. In vitro digestion was based on a combination of the method described by C. Shani-Levi et al. 26 and D. Dupont at al. 27 with several modifications, with the most appreciable change done in the gastric pH profile adjusted to meet that of P. Marteau et al. [START_REF] Marteau | In't Veld[END_REF] who worked on digestion of yoghurt. Simulated gastric fluid (SGF) was freshly prepared by dissolving 2 gr NaCl in 1000ml of double distilled water. Using 2M HCl, the pH was adjusted to 1.2 for pepsin dissolution, or to 5.0 for the initiation of the in vitro digestion model. In practice, the fresh VD 3 -loaded rCM were mixed with SGF (mixing ratio of VD 3 -loaded rCM :SGF 1:1 v/v) and were transferred to 100ml water jacketed reactor (6.1418.250, Metrohm, Switzerland), maintained at 37°C with continuous blending (~250rpm) to mimic the stomach conditions. The initial pH was adjusted to 5, and gastric pepsin (170 U/mg casein, predissolved in SGF pH=1.2) was added to mimic the initial conditions in the stomach. The pH was gradually titrated to 1.8 during 1 hour by a dual auto titration unit (Titrando 902, Metrohm, Switzerland) using 0.5M HCl. The simulated gastric step consisted of five stages; pH was maintained as 5 for 7.5 min, then at pH 4.1, for additional 7.5 min, followed by 7.5 min at pH 3, 7.5 min at pH 2.1, and a final gastric step of 30 min at pH 1.8. The gradual titration was achieved using an auto-titration unit controlled by a dedicated software (Tiamo 2.3 software, Metrohm GA, Switzerland). Samples were then taken after 10, 30 and 60 min. Each sample was quench-frozen using liquid nitrogen to stop the enzyme activity. After 1 hour of simulated gastric step, the digesta underwent 2 hours of an upper intestinal step. Simulated duodenal fluid (SDF) solution was freshly prepared by dissolving 7.0 gr NaCl in 800ml distilled water and the pH was adjusted to 6.5 using 1 M NaOH. Trypsin (34.5 U/mg casein), chymotrypsin (0.4 U/mg casein), pancreatic lipase (85 U/ml), sodium glycodeoxycholate (4mM) and taurocholic acid salt hydrate (4mM) were predissolved in SDF and were added to the digesta to mimic the intestinal conditions. The digesta pH was increased to 6.25 by titration with 1 M sodium bicarbonate solution and maintained constant for 2 hours.

Samples were then taken after 3, 10, 30 and 60 min. Each sample was quench-frozen using liquid nitrogen to stop enzymatic activity.

Protein degradation and VD 3 retention during simulated gastrointestinal digestion:

The degradation of casein as a function of simulated digestion time was studied using sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE). The bands of the protein were stained by coomassie brilliant blue. VD 3 retention as function of digestion time was assessed following the protocol of Kazmi et al. 29 . Briefly, 0.5ml of each sample was transferred into a test tube. A 1.6 ml of mixture of methanol and chloroform (2:1 v/v) was added to each test tube, followed by mixing for 20 seconds. 1 ml of chloroform was then added to each tube, and mixed for 1 min. Samples were then centrifuged for 10 min at 1,500 g at 4°C. From each test tube, 1.2ml of the clear chloroform bottom layer was transferred to an evaporation vial using a glass syringe. The extract was dried from chloroform under a flow of nitrogen gas, reconstituted in 0.75 ml of the reversed phase HPLC (RP-HPLC) mobile phase [methanol: acetonitrile: water (49.5:49.5:1 v:v:v)], and equilibrated for 15 min. VD 3 was quantified by RP-HPLC using HP 1100 system with 4.6 × 250 mm Vydac 201TP C18 column, a diode array detector and an auto-sampler. Twenty microliter samples were run under isocratic flow of 1.3 ml/min. The area under the curve at 267nm (maximal VD absorbance wavelength) was calculated for VD 3 quantification. The ratio between the area under the curve at 228nm and 254 nm was calculated for VD 3 identification.

VD 3 in-vitro bioavailability analysis: uptake by Caco-2 cells

VD 3 uptake experiments: Caco-2 clone TC7 cells were cultured as previously described 30,31 . Cells were seeded and grown in trans-wells for 21 days to obtain confluent and highly differentiated cell monolayers 30,31 . Twelve hours prior to each experiment, the media in apical and basolateral chambers were replaced by serum-free complete medium. At the beginning of the experiment, cell monolayers were washed with 0.5 mL PBS. The apical side of the cell monolayers received either the digested VD 3loaded rCM or the digested free VD 3 diluted in DMEM (1/8), and the basolateral side received serum-free complete medium. After 1 hour of cell exposure with the samples, cells were washed with 0.5 mL ice-cold PBS, then scraped and collected in 0.5 mL PBS. Absorbed VD 3 was estimated as VD 3 present in harvested cells or transported to the basolateral side. All the samples were sealed under nitrogen and stored at -80°C until analysed. VD 3 extraction: VD 3 was extracted from 800 μL aqueous samples using the following method. Distilled water was added to sample volumes below 800 μL to reach a final volume of 800 μL. Ergocalciferol (VD 2 ), routinely used as internal standard, was added to the samples in 800 μL ethanol. The mixture was extracted once with two volumes of hexane. The hexane phase obtained after centrifugation (800 g, 10 min, 4°C) was evaporated to dryness under nitrogen, and the dried residue was dissolved in 200 μL of mobile phase (60% acetonitrile -38% methanol -2% water). A volume of 150 μL was used for HPLC analysis.

VD 3 analysis by HPLC:

The HPLC system comprised a Shimazu separation module (LC-20ADSP HPLC Pumps and SIL-20CHT autosampler, Shimazu, Champs-sur-Marne, France), an SPD-M20A Shimadzu photodiode array detector (PDA, detection at 268 nm, spectral analysis between 190 and 300 nm). VD 3 and VD 2 were separated using a 280×4.6 nm RP C 18 , 8 μm Zorbax column (Interchim, Montluçon, France) and a guard column, as previously described 32 . No VD 3 was detected in basal conditions. Quantification was performed using Chromeleon software (version 6.50 SP4 Build 1000, Dionex) comparing peak area with standard reference curves. All solvents used were HPLC grade from SDS (Peypin, France).

Results and discussion

Extending previous studies on rCM 8,24,25,33 , this study sought better understanding of the technofunctionality and digestive fate of this delivery system and particularly its implications to VD retention and bioavailability.

Light microscopy

Light microscopy images of pure VD 3 in water and VD 3 of the same concentration in the presence of casein (1:1 molar ratio), are presented in Figure 1. VD 3 concentration in the examined samples was 0.1 mg/ml (260 µM). The casein: VD 3 molar ratio was chosen to be similar to the molar ratio we use in rCM.

At this concentration and after an overnight storage at 4°C, it can be seen from Figure 1a that VD 3 forms micron-sized aggregates. However, in the presence of casein, a clear solution is obtained without any visible particles under the light microscope (Figure 1b), indicating submicron co-assembly of the VD 3 and the casein. From this comparison, we can qualitatively conclude that caseins bind VD 3 .

Binding studies by UV absorbance spectra measurements

The UV-absorbance spectrum of casein-only, VD 3 -only, and casein-VD 3 solutions were measured to qualitatively examine whether binding had occurred. The difference between the mathematical summation of UV-absorbance spectrum of the VD 3 -only solution and that of the casein-only solution, and the UV-absorbance spectrum of their solution together is an indication for binding 9 . Observing Figure 2, it can be seen that the mathematical summation curve is lower than that of the mixed casein-VD 3 solution, indicating association between casein and VD 3 . The results indicate that casein and VD 3 associate at this molar ratio (1:1), an assumption supported by the light microscopy images (Figure 1).

Freeze drying and reconstitution

Freeze-drying is widely used for improving chemical and physical stability and extending shelf-life of colloidal systems 34 . During the freezing step there is phase separation into ice and freeze-concentrated solution. The high nanoparticle concentration in this phase may induce irreversible aggregation 34 .

Moreover, adsorption to ice may lead to protein denaturation 35 . To minimize the harmful effect of the freeze-drying process, cryoprotectants are often added. Saccharides are the most popular cryoprotectants in food, as they protect proteins from denaturation, and promote matrix immobilization by raising the glass transition temperature, thereby preventing nanoparticle aggregation 34,36 . After their freeze-drying in the presence of MD, rCM were reconstituted by addition of water and stirring overnight.

We have observed their reconstitution visually and by size distribution measurements by DLS, and the results are presented in Figure 3. The volume-weighted size distributions shows that rCM samples before and after freeze-drying and reconstitution exhibited bimodal distributions (Figure 3). The first subpopulation was comprised of nanoparticles smaller than 100nm and the second was in the range of 100-350 nm. The volume-weighted-size distribution of the rCM before freeze-drying was in line with the previous work performed by our group 8 . It can be seen in Figure 3 that freeze-drying and reconstitution slightly increased the volume weighted percentage of the larger subpopulation, while that of the smaller one decreased accordingly, but shifted to lower sizes. That resulted in a minor visual difference between the samples; it can be seen in the inset of Figure 3 that the reconstituted sample ( 2) is slightly more turbid than the fresh sample (1). However, the process did not significantly affect the average size and it was even slightly lower after reconstitution (average size of 95±2 nm and 89±0.3 nm before and after freeze-drying and reconstitution respectively). About 90% of the VD was retained during the freezedrying and reconstitution process (concentration of 167±8 µg/ml and 151±0.3 µg/ml before and after freeze-drying and reconstitution, respectively). These results indicate that freeze-drying aided by MD, and reconstitution, had only a minor influence on the rCM size distribution, the visual appearance and the concentration of VD.

Simulated gastric and intestinal digestion

In vitro digestion models are often used to screen and interrogate food delivery systems 37 . Various simulated digestion models exist, mimicking the complex physicochemical and physiological conditions in the human gastrointestinal tract 37 . Varying environmental conditions such as pH, ionic strength and enzyme activity induce structural and chemical changes in the delivery system which may trigger bioactive release 37 but might also cause its degradation. To assess casein digestibility and the retention of VD 3 remaining in the lumen of the human gastrointestinal tract, samples of VD 3 loaded rCM and free VD were challenged by a semi-dynamic human digestion model. Results of direct observations, assessment of digestive proteolysis and determination of VD retention are presented in Figures 456. As can be seen in Figure 4, under gastric conditions, VD 3 loaded rCM formed large curd clots, due to the proximity to the isoelectric point of casein (pH=4.6 15 ) and the activity of pepsin 38 . The aggregate size significantly decreased in one hour of simulated gastric digestion, as can be seen in Figure 4. Pepsin activity along with the gradual pH decrease below the isoelectric point, expectedly increased the solubility of the casein as well as its susceptibility to proteolysis. The observed casein clot degradation during the gastric phase was in-line with our SDS-PAGE results shown in Figure 5.

In lane 1 it is possible to see the initial protein composition in the VD 3 loaded rCM suspension. The four casein fractions , α s1 -, α s2 -, β-and κ-caseins were observed in molecular weight range of 19-25 KDa as naturally present in bovine milk 15 . The other bands in lane 1, whose sizes were lower than 19 KDa, were probably residual whey proteins such as β-lactoglobulin and α-lactalbumin (18.3 KDa and 14.2 KDa respectively 15 ) from the sodium caseinate. A few weaker higher-MW bands above 40kDa are also seen, which may be serum albumin, lactoferrin and immunoglobulins, or some casein aggregates which weren't completely dissociated during sample preparation (less likely).

Lanes 2-4 represent the VD 3 loaded rCM under simulated gastric conditions after 10, 30 and 60 min respectively. In these lanes it is possible to witness the degradation of the caseins, as the bands representing the four casein fractions diminished and smaller peptides (2-15 KDa) appeared. The bands whose sizes were between 25 and 37 KDa, were apparently pepsin, but possibly also degradation products of the higher MW bands seen in lane 1. By the end of the simulated gastric phase, the large clot was enzymatically digested and only casein peptides remained (Figure 4 and Figure 5-lane 4). Our results regarding the gastric breakdown of the caseins are in accord with the literature 15,38,[START_REF] Mcsweeney | Advanced dairy chemistry. basic aspects Volume 1A[END_REF] . After 60 min of simulated gastric digestion, the digesta pH was raised to 6.25, and trypsin, chymotrypsin, lipase and bile salts were added to mimic duodenal conditions. Lanes 6-10 (Figure 5) represent 10, 30, 60, 90 and 120 min of simulated upper intestinal digestion. In these lanes it is possible to see bands in the size range of approximately 20-25 KDa which represent the intestinal chymotrypsin, trypsin and lipase that were added in this stage. The addition of duodenal enzymes expectedly caused further casein digestion.

Several bands of casein peptides were observed in lane 6 representing a sample taken after 10 min of simulated duodenal digestion. These bands were not observed in samples taken after 30 min or more of simulated duodenal conditions (lanes 7-10). These SDS-PAGE results are in-line with the visual observations, as the digesta became more transparent with proceeding simulated digestion (Figure 4). In the absence of protein no significant changes were observed over 60 min of simulated gastric digestion, as the vitamin remained dispersed in water and the system remained turbid (Figure 4). During the duodenal phase the system was clarified, most likely due to the addition of bile salts, natural surfactants known to solubilize lipid-soluble compounds and breakdown products, and to form mixed micelles [START_REF] Goodman | [END_REF] . To study VD 3 retention during simulated gastric and duodenal digestion, extraction of VD 3 from the aqueous phase was performed, followed by its quantitative determination by RP-HPLC analysis. Results are presented in Figure 6. It can be seen from Figure 6 that the presence of protein improved VD 3 stability during simulated gastric and intestinal digestion (retention of ~85% vs. ~15% at the end of simulated digestion process in presence and absence of protein respectively). VD 3 is sensitive to the presence of acid conditions and high temperatures, however, VD 3 has better stability under neutral and alkaline conditions [START_REF] Combs | The vitamins[END_REF] . As mentioned above, the initial pH of the simulated gastric phase was 5, which was gradually decreased to 1.8 during one hour. Expectedly, in the absence of protein only 37% of the vitamin survived under these acidic conditions. At the beginning of the simulated duodenal phase the pH was increased to 6.25 and maintained constant for two hours. Under these conditions VD 3 was less susceptible to degradation, and only additional 13% of the free vitamin were degraded in this phase, resulting in 15% of the initial VD 3 remaining after both simulated digestion stages. During simulated digestion of VD 3 loaded rCM, after 10 min of the gastric phase, there was an apparent decrease of the VD 3 residual amount to 40%, however, the apparent residual amount was then increased to 73% and 90% after 30 and 60 min of the gastric phase, respectively. Moreover, the measurement results of these samples were characterized by large error bars. As discussed above, at the beginning of the simulated gastric phase large casein aggregates formed, as seen in Figure 4. It is likely that most of the vitamin was entrapped in these aggregates, thus the samples taken from the bulk during the gastric phase contained less VD 3 and do not accurately represent the overall VD 3 retention status at this stage. The residual amount at the beginning of the duodenal phase, after breakdown and dissolution of the large aggregates, was 90% and it remained relatively constant during this phase. It can be seen in Figure 4 that during the duodenal phase, the casein aggregates were digested and the digesta was clarified, making the sample more homogenous. Thus, it is plausible that the VD 3 was no longer entrapped in casein aggregates but was extracted into mixed micelles and was more uniformly dispersed in the digesta. Therefore, the samples taken at this point better represent the actual VD 3 retention status. Hence, the high values of remaining VD 3 during the duodenal phase, demonstrate the protection conferred to the vitamin by the rCM mainly under simulated gastric digestion conditions. The use of protein nanovehicles to protect VD 3 during simulated digestion was previously examined. Levinson et al. 10 found that 92% of the initial VD 3 concentration survived a two hours simulated gastric step digestion as the VD 3 was complexed within soybean β-conglycinin nanoparticles (relative to 50% retention for VD 3 in phosphate buffer pH 6.8). Lin et al. 11 found that nanocomplexes of corn protein hydrolysate (CPH) with VD 3 with CPH:VD 3 mass ratio of 15:1 and higher conferred protection to the vitamin during simulated digestion (1 hour simulated gastric step followed by 2 hours upper intestinal step), as the remaining vitamin concentration was approximately 95% of its initial concentration 11 . Markman et al. also observed protection of VD by casein and casein-MD Maillard conjugates at pH 2.5 [START_REF] Markman | [END_REF] .

VD 3 in-vitro bioavailability analysis by Caco-2 cells

VD 3 in vitro bioavailability within rCM has been examined using Caco-2 cells. The results are presented in Figure 7. Caco-2 cells were exposed to VD 3 -loaded rCM and to free VD 3 samples taken after simulated digestion including one hour of a gastric step and two hours of upper-intestinal step. Cell exposure time with the digested samples was one hour. This exposure duration was chosen in view of the approximate exposure time of the digesta with a given intestinal epithelial cell, expected to occur in the small intestine considering atypical transit rate 43 , and stomach content/chyme volume 44 (e.g. assuming a typical chyme volume of 1000 ml and intestinal flow rate of 15ml/min, exposure time is ~66 min). The uptake of VD 3 in both formulations is presented in Figure 7A. The VD 3 uptake rates were low for both formulations (approximately 10% and 15% uptake of VD for digested VD 3 -loaded rCM and for digested free VD 3 , respectively). The absorption rate for the digested free VD 3 was slightly, though not significantly higher, compared to that of digested VD 3 -loaded rCM. Taken together with the much higher retention, due to the protective effect of the rCM and the curd-aggregates they form mainly under the acidic gastric conditions, the overall effect of the rCM was that a four-fold higher proportion of the ingested vitamin was bioavailable in vitro, compared to the unencapsulated VD (Figure 7B). These results regarding the high overall effect of rCM are in-line with the high oral bioavailability of VD observed in our preceding clinical studies. It is likely that VD was protected against the harmful acidic environment in vivo by its entrapment within the large protein clots formed at the beginning of the gastric step of digestion (as observed herein), and this can partly explain the high in vivo oral bioavailability obtained.

Conclusions

We studied the performance of VD-loaded rCM after drying and reconstitution, considering VD retention during simulated digestion, and the subsequent in-vitro bioavailability of the vitamin. Our results provided qualitative evidences for binding of VD 3 and caseins, using light microscopy and UV absorbance spectra measurements. Such VD 3 -rCM co-assemblies aided by MD were well reconstituted after freezedrying. The freeze-drying and reconstitution had marginal influence on the size distribution of the VD 3loaded rCM, on their visual appearance and on the VD 3 concentration. Encapsulated VD 3 was also found to be significantly protected during simulated gastric and upper-intestinal conditions, compared with VD 3 dispersed in water. The VD 3 in-vitro bioavailability examined using Caco-2 cells indicated low uptake rates for both formulations, apparently due to limitations of this colon cell-based model (as our clinical studies 8,25,45 showed high bioavailability). However, taken together with the much higher retention of the encapsulated VD 3 during digestion, the overall observed effect of rCM was an improved oral bioavailability of VD 3 compared to the unencapsulated VD 3 . Thus, this work provides yet another support for the great potential of rCM to ensure efficient and more sustainable delivery of sensitive hydrophobic bioactive compounds in food.

Figure captions:

Figure 1 :

 1 Figure 1: VD 3 and casein form submicron co-assemblies. Light microscope images of (a) pure VD in

Figure 2 :

 2 Figure 2: UV absorbance spectra measurements indicating casein-VD 3 binding. The UV-vis absorbance

Figure 3 :

 3 Figure 3: Freeze-drying aided by maltodextrin has only a minor influence on the size distribution and

Figure 4 :

 4 Figure 4: Visual observations of in vitro digestion of free VD3 in water and VD3 loaded rCM show

Figure 5 :

 5 Figure 5: RCM are digested mainly during the gastric step of the digestion process. SDS-PAGE

Figure 6 :

 6 Figure 6: The entrapment in rCM improved VD 3 stability during simulated gastric and intestinal

Figure 7 :

 7 Figure 7: The overall effect of the rCM is a higher proportion of the ingested vitamin being bioavailable

Figure 1 :

 1 Figure 1:

Figure 2 :

 2 Figure 2:

Figure 3 :Figure 4 :

 34 Figure

Figure 5 :Figure 6 :

 56 Figure

Figure 7 :

 7 Figure 7:

Acknowledgements

The authors would like to thank Mrs. Maya Bar-Zeev and Dr. Carmit Shani-Levi for their technical assistance.


Olivier Beaumont

Brett A Becker

Ashley Deflumere

Lionel Eyraud-Dubois

Thomas Lambert

Alexey Lastovetsky

Recent Advances in Matrix Partitioning for Parallel Computing on Heterogeneous Platforms

come   L'archive ouverte pluridisciplinaire

B is computed. Let us assume that at each step, processor P computes s blocks of matrix C, whose projections along the different axes have respective sizes h and w. Then, the volume of computation for P is proportional to s and the volume of communication required for P to complete its computation is proportional to h + w, since it needs to receive the corresponding blocks of both matrices A and B. In order to balance the computing load, each processor should compute a number of blocks of C proportional to its relative speed. In turn, the overall volume of communication is proportional to the sum of the projections of the areas owned by the different processors along the axes. Therefore, in order to minimize the processing time while also minimizing the overall volume of communication, the optimization problem is amenable to the problem of partitioning a square into a set of zones of prescribed area (in order to balance the load) such that the sum of the projections along the two axes is minimized (in order to minimize the communications).

This paper contrasts results of two approaches to this optimization problem when partitioning among small numbers of heterogeneous resources. The first approach is that of Becker and Lastovetsky who proposed that non-rectangular partitionings can be optimal for partitioning between two heterogeneous resources [START_REF] Becker | Matrix multiplication on two interconnected processors[END_REF], and proven to be optimal among all possible partitionings for two resources (for certain power ratios) by DeFlumere, Lastovetsky and Becker [START_REF] Deflumere | Partitioning for parallel matrix-matrix multiplication with heterogeneous processors: The optimal solution[END_REF]. This approach was extended to three resources [START_REF] Becker | Towards data partitioning for parallel computing on three interconnected clusters[END_REF], and proven (again for certain power ratios) in [START_REF] Deflumere | Optimal partitioning for parallel matrix computation on a small number of abstract heterogeneous processors[END_REF]. The second approach is based on recursive approximation algorithms. Recently, Beaumont, Eyraud-Dubois and Lambert improved the best-known approximation ratio to 2 √ 3 ≈ 1.15 [START_REF] Beaumont | A new approximation algorithm for matrix partitioning in presence of strongly heterogeneous processors[END_REF] by mixing the work on optimal partitionings for small numbers of partitions and approximation algorithms for rectangle partitioning proposed by Nagamochi et al. [START_REF] Nagamochi | An approximation algorithm for dissecting a rectangle into rectangles with specified areas[END_REF] and Fügenschuh [START_REF] Fügenschuh | Exact and approximation algorithms for a soft rectangle packing problem[END_REF].

II. GENERAL CONTEXT

A. Related Works

This optimization problem was first introduced by Lastovetsky and Kalinov in [START_REF] Kalinov | Heterogeneous distribution of computations solving linear algebra problems on networks of heterogeneous computers[END_REF]. In [START_REF] Beaumont | Partitioning a square into rectangles: NP-completeness and approximation algorithms[END_REF], it was proven that the problem is NP-Complete, and a first approximation algorithm with a bounded ratio of 1.75 was proposed. The assumption of rectangular partitions was relaxed in [START_REF] Becker | Matrix multiplication on two interconnected processors[END_REF] for partitioning between two resources. This non-rectangular partitioning has perfect load balance and a lower overall volume of communication compared to all rectangular partitionings provided the computational lpower ratio between the two resources exceeds 3 : 1. This was extended to three resources in [START_REF] Becker | Towards data partitioning for parallel computing on three interconnected clusters[END_REF] where a non-rectangular partitioning based on the two partition case was proposed which again has perfect load balancing, and a lower overall volume of communication, again subject to certain power ratio and topology requirements. It is important to note that these partitioning algorithms are not modifications of homogeneous ones; they are designed for heterogeneous scenarios [START_REF] Becker | High-level data partitioning for parallel computing on heterogeneous hierarchical computational platforms[END_REF]. These algorithms have also been shown to benefit real-world applications [START_REF] Becker | Max-plus algebra and discrete event simulation on parallel hierarchical heterogeneous platforms[END_REF]. It has also been demonstrated that optimizing the partitioning for non-square matrices is not achievable with straight-forward scaling of the square cases, at least for specific applications such as matrix multiplication [START_REF] Becker | High-level data partitioning for parallel computing on heterogeneous hierarchical computational platforms[END_REF].

In [START_REF] Deflumere | Partitioning for parallel matrix-matrix multiplication with heterogeneous processors: The optimal solution[END_REF] it was proven that in the case of two resources, the optimal partitioning for any power ratio is of one of two forms, one rectangular and one non-rectangular. A method called the Push technique was first proposed in this work in order to generate a small number of potentially optimal shapes, which can then be analyzed to find the optimal one. The Push technique was extended and used in the study of three resources in [START_REF] Deflumere | Searching for the optimal data partitioning shape for parallel matrix matrix multiplication on 3 heterogeneous processors[END_REF] in order to identify potential optimal shapes. Six potential shapes were identified in this work. While the optimality of the three resource case has never been mathematically proven, extensive simulations conducted in [START_REF]Optimal data partitioning shape for matrix multiplication on three fully connected heterogeneous processors[END_REF] and [START_REF] Deflumere | Optimal partitioning for parallel matrix computation on a small number of abstract heterogeneous processors[END_REF] failed to find a shape that would outperform these six shapes.

Independently, recursive partitioning algorithms for this optimization problem have been recently proposed. In these, the set of processors is split into two parts at each step, converging on an approximation of the optimal partitioning. Sophisticated proof techniques enabled Nagamochi and Abe [START_REF] Nagamochi | An approximation algorithm for dissecting a rectangle into rectangles with specified areas[END_REF] to improve the approximation ratio down to 1.25. Recently, Fügenschuh et al. [START_REF] Fügenschuh | Exact and approximation algorithms for a soft rectangle packing problem[END_REF] improved this result to 1.15, but under the assumption that if we consider processors in decreasing order of their processing speeds, there is no abrupt change in the performance between any two successive processors. Unfortunately, such an abrupt decrease typically happens when considering more heterogeneous nodes such as those consisting of CPUs and GPUs, therefore restricting Fügenschuh's algorithm to the case of loosely heterogeneous platforms.

Beaumont, Eyraud-Dubois and Lambert [START_REF] Beaumont | A new approximation algorithm for matrix partitioning in presence of strongly heterogeneous processors[END_REF] sought to keep the best of both worlds by adapting the idea of non-rectangular partitionings proposed by Becker and Lastovetsky and extending it to an arbitrary number of processors by adapting the recursive partitioning algorithm proposed by Nagamochi, which facilitates approximation-ratio proofs. These two ingredients led to an improvement of the approximation ratio down to 2 √ 3 ≈ 1.15. This non-rectangular recursive partitioning algorithm (NRRP), does not require any specific assumption on the relative speed of resources and is therefore applicable to computational power ratios comparable to those of nodes consisting of both regular cores and accelerators.

This partitioning problem can be used as a building block for many dense linear algebra kernels. For instance, it has been extended to LU factorization and other dense linear algebra kernels in [START_REF] Beaumont | Static LU decomposition on heterogeneous platforms[END_REF], [START_REF] Beaumont | A proposal for a heterogeneous cluster ScaLAPACK (dense linear solvers)[END_REF]. In this case, a block cyclic principle is combined with the initial partitioning in order to obtain 2D cyclic ScaLAPACK solutions [START_REF] Blackford | ScaLAPACK users' guide[END_REF], where the load is balanced throughout the whole computation. These partitionings have also been adapted to distributed hierarchical and highly heterogeneous platforms in [START_REF] Clarke | Hierarchical partitioning algorithm for scientific computing on highly heterogeneous CPU + GPU clusters[END_REF], where the partitioning is applied at two levels (intra-node and internode), based on sophisticated performance models. The same partitioning has also been extended to finite-difference timedomain (FDTD) methods to obtain numerical solutions of Maxwell's equations in [START_REF] Shams | On optimization of finite-difference timedomain (FDTD) computation on heterogeneous and GPU clusters[END_REF]. The extension to more dynamic settings has also been considered in [START_REF] Mohamed | DDOps: dual-direction operations for load balancing on non-dedicated heterogeneous distributed systems[END_REF]. In this case, the partitioning problem can be used in order to provide an initial static partitioning algorithm that can be modified in order to dynamically maintain load balancing. Recently, in order to cope with resource heterogeneity and the difficulty in building optimal schedules, the use of dynamic runtime schedulers has been proposed, such as StarPU [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF], StarSs [START_REF] Planas | Hierarchical task-based programming with StarSs[END_REF], QUARK [START_REF] Yarkhan | QUARK users' guide: Queueing and runtime for kernels[END_REF] or PaRSEC [START_REF] Bosilca | PaRSEC: Exploiting heterogeneity to enhance scalability[END_REF]. At runtime, the scheduler takes the scheduling and allocation decisions based on the set of ready tasks (tasks for which data and control dependences have been solved), on the availability of the resources (estimated using expecting processing and communication times), and on the actual location of input data. The comparison between static scheduling strategies (such as the one proposed in this paper) and runtime scheduling strategies has been recently considered in [START_REF] Beaumont | Comparison of static and dynamic resource allocation strategies for matrix multiplication[END_REF], where the analysis of the behavior of static, dynamic, and hybrid strategies highlights the benefits of introducing more static knowledge and allocation decisions in runtime libraries.

The results of all of these papers are based on the partitioning problem considered in this paper and can therefore directly benefit from an improvement in the performance and approximation ratio.

B. Notations and Problem Statement

In this section, we define the notations that will be used in the rest of this paper and we present the formal version of the optimization problem that corresponds to enforcing a perfect load balancing while minimizing the amount of communications.

The problem comes in several flavors depending on the constraints imposed to the partition. In this paper, in order to have a solution that is independent of the matrix size (but only dependent on the relative speeds of the resources), we will consider the continuous version of the problem, where the goal is to split the unit square [0, 1] × [0, 1]. In general, such a solution has to be rounded for a given matrix size. How to perform these roundings optimally has been considered in [START_REF] Beaumont | Comparison of static and dynamic resource allocation strategies for matrix multiplication[END_REF] and we rely on these results. Another possible restriction is related to the shape of the area Z l allocated to processor P l . In many of the above mentioned works, it has been assumed that Z l had to be shaped as a rectangle [START_REF] Beaumont | Partitioning a square into rectangles: NP-completeness and approximation algorithms[END_REF], [START_REF] Nagamochi | An approximation algorithm for dissecting a rectangle into rectangles with specified areas[END_REF], [START_REF] Fügenschuh | Exact and approximation algorithms for a soft rectangle packing problem[END_REF] because this was considered as a requirement for a simple and efficient implementation. Recent task-based implementations allow this constraint to be relaxed, and it is thus possible to search for general (non-rectangular) partitionings. Such partitionings have been proposed in [START_REF]Optimal data partitioning shape for matrix multiplication on three fully connected heterogeneous processors[END_REF], [START_REF] Becker | Towards data partitioning for parallel computing on three interconnected clusters[END_REF], and it has been proven that it is possible to derive significantly better partitionings when removing the rectangular constraint. Therefore, in this paper, we do not impose that the areas allocated to the processors be shaped as rectangles.

We thus consider the unit square S = [0, 1] × [0, 1]. Let Z denote a zone (a subset of S) included in the unit square. We denote by s(Z) its area (formally Z dxdy) and by R(Z) its covering rectangle, i.e. the Cartesian product of the projections of Z along both dimensions.

If R(Z) = [x 1 , x 2 ] × [y 1 , y 2 ],
then we define the height of Z by h(Z) = x 2 -x 1 and the width of Z by w(Z) = y 2 -y 1 . Finally, let us define p(Z) = h(Z) + w(Z), the half-perimeter of R(Z) and ρ(Z) = max(h(Z), w(Z)) min(h(Z), w(Z)) its aspect ratio.

We consider the following problem :

Problem 1 (PERI-SUM). Given a set of n positive rational numbers {s 1 , . . . , s n } such that s k = 1, and the square S = [0, 1] × [0, 1], find for each s k a zone Z k ⊆ S such that the area of Z k is s k , Z k = S, and such that p(Z k ) is minimized.

In the following, we denote p(Z k ) as c(Z 1 , . . . , Z n ) and its optimal value as c opt . A lower bound has been proposed by Ballard et al. in [START_REF] Ballard | Minimizing communication in linear algebra[END_REF], that comes from an application of the Loomis-Whitney inequality. This lower bound simply states that the perimeter of a zone Z k of given area s(Z k ) is minimal when the zone is shaped as a square.

c(Z k ) ≥ 2 s(Z k ) (1)
This provides a lower bound on the total cost for any instance, c opt ≥ 2 k √ s k . This bound is reached only when all zones are squares, which may not be feasible depending on the instance (consider for instance the case of two identical zones of area 1 /2), and hence this lower bound is in general optimistic.

C. Complexity Results

A variant of the decision problem associated to this optimization problem has been proved to be NP-Complete. More specifically, it has been proven in [START_REF] Beaumont | Partitioning a square into rectangles: NP-completeness and approximation algorithms[END_REF] that deciding whether S = [0, 1] × [0, 1] can be partitioned into squares of respective areas {s 1 , . . . , s p } is NP-Complete.

Since squares are the best possible shapes for the zones irrespective to any additional constraint, this complexity result automatically translates to PERI-SUM, which is thus an NP-Complete problem. For this reason, the rest of the paper is devoted to finding algorithms in two directions: either exact algorithms for a small number of processors, or approximation algorithms for the general case.

D. Theoretical comparison with block-cyclic approaches

Distributed matrix multiplication implementations often use classic partitioning schemes, such as 2D cyclic or blockcyclic decompositions. In such decompositions, processors are arranged into a p × q grid: all processors on the same row receive the same rows of matrix A, and all processors on the same column receive the same columns of matrix B. These decompositions can only be used in the (very) special case of homogeneous processors, since they are inherently designed to provide each processor with the same amount of computation. Designing block-cyclic decompositions that would best adapt to heterogeneous computing speeds is itself a difficult problem. Since we focus on heterogeneous settings, we can not compare our results directly to block-cyclic decompositions.

For homogeneous processors, since each element of matrix A is sent to p processors, and each element of matrix B is sent to q processors, the total communication cost of a blockcyclic decomposition is p+q. When p = q (= √ n), the cost is 2 √ n, which is exactly equal to the lower bound 2

n i=1 1 √
n . However, if p and q are very different, the ratio to the lower bound can be very large. The worst case is reached when one of p or q is equal to 1, where the communication cost is n + 1, which yields an unbounded ratio with respect to the lower bound.

III. OPTIMAL PARTITIONINGS

Traditionally many applications such as matrix multiplication (which is generalizable to many other applications of interest) are performed in parallel by dividing matrices into rectangular partitionings. We challenge this approach. In this section we study the optimality of partition shapes for a small number of partitions. Traditionally the problem of finding the optimal shape is reduced to PERI-SUM, the problem of finding the optimal partitioning that minimizes the sum of half perimeters (SHP), which, as discussed in Section II, is equivalent to the projections of the areas owned by the different processors along the axes. This is the simplest mathematical formulation of this problem. Our motivation is that when we balance the load (which is normally trivial) and minimize the SHP, we minimize the amount of data moved between partitions, and thus the required communication for a given application.

For small numbers of heterogeneous processors, potentially optimal partition shapes can be determined using the Push Technique [START_REF] Deflumere | Partitioning for parallel matrix-matrix multiplication with heterogeneous processors: The optimal solution[END_REF], [START_REF] Deflumere | Searching for the optimal data partitioning shape for parallel matrix matrix multiplication on 3 heterogeneous processors[END_REF].

In the following, we denote by P the fastest processor, by Q the second fastest, and when applicable, by R the slowest one (s P ≥ s Q ≥ s R ).

A. Two Heterogeneous Processors

Two partition shapes shown in Figure 1 are identified as the optimal two processor partitioning in [START_REF] Deflumere | Partitioning for parallel matrix-matrix multiplication with heterogeneous processors: The optimal solution[END_REF].

Q P Square Corner Q P Straight Line
Figure 1: The optimal partitions with two processors.

Theorem 1. For two heterogeneous processors, the data partitioning with the optimal total volume of communication will be the Square Corner when the ratio of computational power between the two processors is greater than three.

Proof. The fact that the two partitions of Figure 1 are dominant is a corollary of Theorem 2, proven below. We can compute the cost of each partition:

• Square Corner. The cost for P is 2, and the cost for

Q is 2 √ s Q , for a total cost of c SC = 2 + 2 √ s Q .
• Straight Line. The cost for P is 1 + s P , and the cost for

Q is 1 + s Q , with s P + s Q = 1, for a total cost of c SL = 3.
Alternatively, we can see that one dimension is sent to both processors P and Q, while the other is shared among them, hence the cost of 3.

It is clear that c SC < c SL if and only if √ s Q < 1 2 , which is equivalent to s Q < 1 4 .

B. Three Heterogeneous Processors

In [START_REF] Deflumere | Searching for the optimal data partitioning shape for parallel matrix matrix multiplication on 3 heterogeneous processors[END_REF], the Push technique was used to propose six candidate partition shapes for consideration to be the optimal three processor partitioning of a square shown in Figure 2. The optimality of these shapes has never been proven however. What was proven is that the Rectangle Corner and the Straight Line both have a higher communication cost than the Block Rectangle [START_REF]Optimal data partitioning shape for matrix multiplication on three fully connected heterogeneous processors[END_REF]. It has also become evident that the L Rectangle and the Block Rectangle are topologically equivalent. In this paper, we prove that the four shapes shown in Figure 3 are sufficient to optimally partition an arbitrary rectangle. We also prove that three of these shapes -the Square Corner, the Square Rectangle, and the Block Rectangle, are sufficient to optimally partition a square.

Theorem 2. The optimal data partition shape to split a rectangle between three processors is always one of the four shapes of Figure 3: the Square Corner, the Square Rectangle, the Block Rectangle, or the Straight Line.

Proof. Without loss of generality, we assume that the rectangle to partition has dimension a × b, with a ≤ b, and a is the size of the vertical axis. We first compute the cost of each solution (note that the first two might not be feasible):

R Q P Square Corner a b P Q R Square Rectangle P Q R Block Rectangle P R Q Straight Line
• Square Corner c SC = a + b + 2( √ s Q + √ s R ) • Square Rectangle c SR = 2a + b + 2 √ s R • Block Rectangle c BR = 2a + 2b -sP a • Straight Line c SL = 3a + b
To prove Theorem 2, we consider any valid partition S, and we prove that the cost C of S is at least as much as the cost of one of the proposed shapes. The proof is split into 5 different cases, depending on the number of resources that hold a complete row or column in S. In the following, we say that a processor holds a complete row (resp. column) if the projection of its allocated area on the horizontal (resp. vertical) axis is equal to the complete side of the rectangle.

(i) No processor holds a complete side in S.

In this case, each data is sent to at least 2 resources, so that the overall size of the projections is at least 2a + 2b.

Since the Straight Line position in Figure 3 is always feasible and has cost 3a + b ≤ 2a + 2b, S is dominated. (ii) Exactly one processor T holds a complete side in S, but not both sides.

Let us first assume that processor T holds a complete column, and let us denote by x the length of the part of the horizontal axis on which only the area allocated to T is projected (i.e., the amount of horizontal data that is sent only to T , see Figure 4). Then, the other b -x columns are sent to at least 2 resources (otherwise another resource would hold a complete side), and each row is sent to at least 2 resources as well (T and at least another one) so that the cost of S is at least

2a + x + 2 * (b -x) = 2b + 2a -x. By construction x ≤ s T a
, and the cost of the Block Rectangle solution is 2b + 2a -s P a , with s P ≥ s T . Hence S is dominated. Similarly, if T holds a complete row, the cost of S is at least 2a+2b-x, with x ≤ s T b , and is thus also dominated by Block Rectangle since a ≤ b.

(iii) Exactly one processor T holds a complete row and a complete column in S.

In this case, T contributes for at least a + b to the cost of S, and therefore

C ≥ a + b + 2 √ s U + 2 √ s V , where
U and V are the two other processors. For example, let us assume that they both hold a complete column. Since the remaining resource V does not hold a complete column, each data on the horizontal axis is sent at least to either T or U . Then, the contribution of T and

• If √ s U + √ s V > a, then C > 3a + b and S is again dominated. • If √ s U + √ s V < a,
U to C is at least 2a+b and therefore C ≥ 2a+b+2 √ s V .
If both resources hold a complete row instead, we obtain similarly that

C ≥ 2b + a + 2 √ s V ≥ 2a + b + 2 √ s V . • If √ s V > a, then C > 4a + b. Since the Straight
Line solution has cost 3a + b, this implies that S is dominated (by Straight Line solution). It is interesting to note that the Straight Line solution is never optimal when partitioning a square (i.e., when a = b), since it is always dominated by Block Rectangle. Of the remaining three, the first two are non-rectangular, but nonrectangular in different ways which make each suited to a different type of heterogeneous distribution of computational power. In general, these results show that, for a square partitioning (see Figure 11 in Section V):

• If √ s V < a, then solution Square Rectangle from
• Square Corner is the optimal shape for heterogeneous systems with a single fast processor, and two relatively slow processors. • Square Rectangle is the optimal shape for heterogeneous systems with two fast processors, and a single relatively slow processor. • Block Rectangle is the optimal shape for heterogeneous systems with a fast, medium, and slow processor, as well as relatively homogeneous systems.

IV. NON-RECTANGULAR RECURSIVE PARTITIONING (NRRP) Extending the Push Technique to an arbitrary number of partitions is not trivial, and the challenge increases significantly for each additional partition (the three partition case was much more challenging than the two partition case). This creates a demand for approximation algorithms for arbitrary numbers of partitions. We are inspired by the recursive results of Nagamochi and Abe [START_REF] Nagamochi | An approximation algorithm for dissecting a rectangle into rectangles with specified areas[END_REF] whose recursive algorithm reduced the approximation ratio to 1.25 without significantly restricting the problem. We develop a recursive algorithm that guarantees not to be more than 2 √ 3 from optimal. Our motivation is to be able to compute fast a partitioning that is known to be within a small factor from the optimal for any arbitrary number of partitions.

This section is devoted to approximation algorithms for PERI-SUM. We present several increasingly complex algorithms, to finish with NRRP (Non-Rectangular Recursive Partitioning), whose approximation ratio is below 2 √ 3 .

A. Rectangular Recursive Partitioning

Algorithm 1: RRP(R, {s 1 , . . . , s n })

Input: A rectangle R, a set of positive values {s 1 , . . . , s n } such that s i = s(R) and

s 1 ≤ s 2 ≤ . . . ≤ s n Output: For each 1 ≤ i ≤ n, a rectangle R i such that s(R i ) = s i and R i = R if n = 1 then return R else ρ = ρ(R) ; k = the smallest k such that k i=1 s i ≥ s 3 ; s = k i=1 s i ; Cut R in R 1 ,
R 2 according to its greatest dimension and with s(R 1 ) = s ; return RRP(R 1 , {s 1 , . . . , s k }) + RRP(R 2 , {s k+1 , . . . , s n })

The first algorithm is called RRP (Rectangular Recursive Partitioning, see Algorithm 1) and is inspired by the seminal work of Nagamochi et al. [START_REF] Nagamochi | An approximation algorithm for dissecting a rectangle into rectangles with specified areas[END_REF]. The basic idea is to use the divide and conquer paradigm by splitting the current rectangle in two sub-rectangles at each step of the algorithm, and assign to each part a sublist of the list of processors. The cut is made so as to have, if possible, at least a third of the total area on each part. For this purpose, the areas are sorted in increasing order and aggregated one by one until one third of the total is reached. Lemma 3 shows that if at least one item remains, then this method ensures that both parts contain at least one third of the total area. Lemma 3. Let {s 1 , . . . , s n } be a set of positive values such that s 1 ≤ . . . ≤ s n and k be the smallest k such that

k i=1 s i ≥ s 3 . Then, if k < n, n i=k+1 s i ≥ s 3 .
Proof. By definition of k,

k-1 i=1 s i < s 3 .
Let us assume that n k+1 s i < s 3 for the search of a contradiction. In this case we obtain that s k ≥ s 3 . As s k+1 ≤ n i=k+1 s i < s 3 , we have s k > s k+1 which is a contradiction with the hypothesis

s 1 ≤ . . . ≤ s n . Otherwise, if k = n, this means that n-1 i=1 s i < s
3 , and hence s n is significantly greater than the other values. This is actually the pathological case on which the following algorithms improve. If we disregard this case for a moment, then the algorithm is quite effective and simply using Lemma 4 leads directly to an approximation ratio of 2 √ 3 , as Fügenschuh et al point out [START_REF] Fügenschuh | Exact and approximation algorithms for a soft rectangle packing problem[END_REF]. Note that the proof of Lemma 4 relies on the lower bound from (1) and does not compare directly to the true optimal value. Lemma 4. Let R be a rectangle. Then:

p(R) 2 s(R) = 1 + ρ(R) 2 ρ(R) .
In particular, if ρ(R) ≤ 3 then:

p(R) 2 s(R) ≤ 2 √ 3 .
Proof. Let us assume without loss of generality that

h = h(R) ≤ w = w(R) and denote ρ = ρ(R) = w h . Then p(R) = h + w = h(1 + ρ). In addition s = s(R) = hw = ρh 2 . Therefore w + h 2 √ s = h(1 + ρ) 2 ρh 2 = 1 + ρ 2 √ ρ . Furthermore, one can prove that x → 1+x 2 √
x is an increasing function on [START_REF] Becker | Matrix multiplication on two interconnected processors[END_REF][START_REF] Becker | Towards data partitioning for parallel computing on three interconnected clusters[END_REF] and then, for ρ ≤ 3:

w + h 2 √ s ≤ 1 + 3 2 √ 3 = 2 √ 3 .
However, the pathological cases imply that the actual approximation ratio of RRP is at least 3 2 . Indeed, let us consider an instance with two values and 1 -. The optimal partitioning is depicted on the left of Figure 5, the communication cost is 2 √ for the zone of area and 2 for the other zone, which yields a total of 2(1 + √ ). The partitioning returned by RRP is shown on the right of Figure 5. In this case, the communication cost is 1 + for the smaller zone and 2for the larger one, which gives a total communication cost of 3. It results that Cost(RRP

) Copt = 3 2(1+ ) -→ →0 3
2 . In their work, Nagamochi et al. prove that RRP is actually a 5 4 -approximation algorithm for the variant of the problem where all the zones have to be rectangles (and indeed in this case the partitioning on the right of Figure 5 is the only possible partitioning and therefore also the optimal one).

1 - √ 1 √ 1 √ 1 - 1 - 1 rien

B. Simple Non-Rectangular Recursive Partitioning

As a first step to improve RRP, we use the insight from Section III: whenever a value is significantly bigger than the others, it is best to avoid splitting into two rectangles. We thus slightly adapt the algorithm and obtain two subroutines.

Guillotine: The first one is the Guillotine routine, depicted in Figure 6 (left). It is the main ingredient of RRP. Given a composed zone R and a rational number α ∈ [0, 1], Guillotine(R, α) splits R along the largest dimension into two rectangles of respective areas αs(R) and

(1 -α)s(R). R 1 R 2 R R 1 Z 2 R

Square:

The second routine is the Square routine, depicted in Figure 6 (right). Given a rectangle R and a rational number α ∈ [0, 1], Square(R, α) returns a square R 1 of area αs(R) and a zone Z 2 which corresponds to the initial rectangle R punched by square R 1 . The covering rectangle of Z 2 is R and Z 2 will always be used to host a simple zone.

With these two subroutines we propose an improved algorithm: SNRRP (Simple Non-Rectangular Recursive Partitioning), given in Algorithm 2. The basic idea is the following: according to an appropriate condition, the current rectangle is either split into two well-shaped rectangles (lines 7-9), or into a square and its complement, the latter being assigned to a single processor (lines 10-12).

Algorithm 2: SNRRP(R, {s 1 , . . . , s n })

Input: A rectangle R, a set of values {s 1 , . . . , s n } such that

s i = s(R) and s 1 ≤ s 2 ≤ . . . ≤ s n Output: For each 1 ≤ i ≤ n, a zone Z i such that s(Z i ) = s i and Z i = R 1 if n = 1 then 2 return R 3 else 4 ρ = ρ(R) ; 5 k = the smallest k such that k i=1 s i ≥ s 3ρ ; 6 s = k i=1 s i ; 7 if k < n then 8 R 1 , R 2 = Guillotine(R, s /s) ; 9 return SNRRP(R 1 , {s 1 , . . . , s k }) + SNRRP(R 2 , {s k+1 , . . . , s n }) 10 else 11 R 1 , Z 2 = Square(R, (s -s n )/s) ; 12 SNRRP(R 1 , {s 1 , . . . , s n-1 }) + Z 2
This algorithm can be proven to be a 1.22) and thus be seen as a significant improvement of RRP. The proof relies on one major invariant (that also ensures the correctness of the algorithm): each rectangle on which the function is called has an aspect ratio below 3, as is steted in Theorem 5.

Theorem 5 (Correctness).

When executing SNRRP(R, {s 1 , . . . , s n }) with ρ(R) ≤ 3, all the recursive calls to SNRRP(R , {s 1 , . . . , s k }) are performed on a rectangle area R such that ρ(R ) ≤ 3.

Proof. Recursive calls to SNRRP are made on two kinds of rectangles, those produced by Guillotine (line 8) and the square produced by Square (line 11). In the second case the aspect ratio is trivially below 3. In the first case, let s and s be as described in Algorithm 2, α = s /s and let R 1 , R 2 = Guillotine(R, α). We have to prove that ρ(R 1 ) ≤ 3 and ρ(R 2 ) ≤ 3. By hypothesis we know that α ≥ 1 3ρ(R) .

Let us assume without loss of generality that

h = h(R) ≤ w = w(R) and denote ρ = ρ(R) = w h . By definition h(R 1 ) = h(R) and w(R 1 ) = αw. Therefore ρ(R 1 ) = max( αw h , h αw ). We have αw h ≤ w h = ρ ≤ 3. In addition, since α ≥ 1 3ρ , h αw ≤ 3ρh w = 3. Hence, if α ≥ 1 3ρ(R) then ρ(R 1 ) ≤ 3.
To prove that ρ(R 2 ) ≤ 3 we can slightly adapt Lemma 3 and show that 1 -α ≥ 1 3ρ(R) . Then, similar arguments prove that ρ(R 2 ) ≤ 3 and yield the proof of the Theorem.

To achieve the approximation proof, let us first note that the zones returned by the algorithm are produced at line 2 and 12. In the first case, the zone is a rectangle with an aspect ratio below 3 (Theorem 5) and Lemma 4 proves that its half-perimeter is below 2 √ 3 times the lower-bound (and

2 √ 3 ≤ 3
2 ). In the case of line 12, the resulting zone is no longer rectangular and an additional lemma is required.

Lemma 6. Let R be a rectangle such that ρ(R) ≤ 3, α ∈ [0, 1] and R , Z = Square(R, α). If α ≤ 1 3ρ(R) then p(Z) 2 √ s(Z) ≤ 3 2 .
Proof. Let us denote ρ(Z) = ρ and s(Z) = s and let us suppose that h = h(R) ≤ w = w(R) without loss of generality. Note that the covering rectangle of Z is R and therefore w = w(Z) and h = h(Z). In this case w = ρh, which implies p(Z)

= (1 + ρ)h. By definition of Square, s = (1 -α)s(R) and s(R) = hw = ρh 2 . Therefore s = (1 -α)ρh 2 ≥ (1 -1 3ρ(Z) )ρh 2 and hence p(Z) 2 S(Z) ≤ (1 + ρ)h 2 (1 -1/3ρ)ρh 2 ≤ √ 3(1 + ρ) 2 (3ρ -1) .
One can prove that the function x -→ √ 3(1+x) 2 √ 3x-1 reaches its maximum on [START_REF] Becker | Matrix multiplication on two interconnected processors[END_REF][START_REF] Becker | Towards data partitioning for parallel computing on three interconnected clusters[END_REF] for x = 1 and x = 3 and this maximum is 3 2 . Thus, p(Z)

2 √ S(Z) ≤ 3 2 .
Lemma 6 covers the second case of zone creation. Indeed, the hypothesis of the lemma on the aspect ratio of R is ensured by Theorem 5 and the hypothesis on α comes from the fact that, in line 11 of Algorithm 2, s -

s n = n-1 i=1 s i = k-1 i=1 s i < s 3ρ(R) by definition of k in line 5.
All of the above proves that if {Z 1 , . . . , Z n } = SN RRP (R, {s 1 , . . . , s n }), then for all i ∈ [1, n], p(Zi)

2 √ si ≤ 3 2 . Therefore, n i=1 p(Zi) n i=1 2 √ si ≤ 3 2 . Since n i=1 2
√ s i is the lower bound from (1), this achieves the proof Theorem 7.

Theorem 7. SNRRP is a 3 2 -approximation for PERI-SUM. Note that there are cases where Algorithm 2 returns a partition such that the ratio between the sum of perimeters and the lower bound based on Equation 1 is indeed 3 2 . Let us define for 1 < k ≤ n, s k = 2/3 n-k+1 and s 1 = 1/3 n-1 . One can notice that s k = 2 3 k i=1 s i . Hence each step of SN RRP ([0, 1] 2 , {s 1 , . . . , s n }) corresponds to the case of line 11-12 and therefore the Square routine is called. In addition, this corresponds the extremal position of the case of line 11-12 and we can notice in the proof of Lemma 6 that in this case the bound is tight. Hence, if Z 1 , . . . , Z n = N RRP ([0, 1] 2 , {s 1 , . . . , s n }), p(Z 1 ) = 2 s(Z 1 ) and for

k > 1, p(Z k ) = 2 3 2 s(Z k ). Thus, lim n→∞ n i=1 p(Z i ) 2 n i=1 s(Z i ) = 3 2
An illustration of this case is depicted in Figure 7.

C. Non-Rectangular Recursive Partitioning

Further improvements over SNRRP and its approximation ratio can be obtained by changing the invariant. In both RRP and SNRRP, the aim is to obtain a rectangle with an aspect ratio below 3. This can be generalized by aiming for a rectangle with an aspect ratio below µ where µ is a parameter whose value will be determined later. The objective is to use values of µ < 3, so that the zones produced are as close to squares as possible; the algorithm NRRP that we propose uses µ = 5 2 . With values of µ different from 3, the proofs change at two places.

1) Aspect ratio guarantees: The first aspect of the proof of SNRRPis Lemma 3 which ensures the invariant that all recursive calls are made on rectangles with aspect ratio below 3. This lemma can be generalized for µ ≥ 3 and can be rewritten as Lemma 8.

Lemma 8. For µ ≥ 3, let {s 1 , . . . , s n } be a set of positive values which sum to s such that s 1 ≤ . . . ≤ s n , and let k be the smallest k such that

k i=1 s i ≥ s µ . Then, if k < n, n i=k+1 s i ≥ s µ .
However, in the case µ < 3, the lemma has to be significantly weakened, as expressed in Lemma 9.

Lemma 9. Let p > 1 be an integer and let µ be in [ 2p+1 p , 2p-1 p-1 [. Let {s 1 , . . . , s n } be a set of positive values such that s 1 ≤ . . . ≤ s n and k be the smallest k such that

k i=1 s i ≥ s µ . Then, if k ≤ n -p, n i=k+1 s i ≥ s µ . Proof. By definition of k, k-1 i=1 s i < s µ .
Let us assume n i=k+1 s i < s µ for the search of a contradiction. In this case we obtain that s k ≥ µ-2 µ s. In addition

n i=k+1 s i ≥ n i=k+1 s k+1 = (n -k)s k+1 . Therefore, s k+1 < s (n-k)µ . From the assumption k ≤ n -p we get that s k+1 < s pµ . Furthermore, µ -2 ≥ 2p+1 p -2 ≥ 1 p . Hence s k ≥ s pµ > s k+1 , what is a contradiction with s 1 ≤ . . . ≤ s n .
Both lemmas provide a sufficient condition for the possibility to split the list of values such that each part sums to a fraction of at least 1 µ of the total, allowing to satisfy the invariant when performing a recursive call on each sublist. When the condition is not met, the implication in the case µ ≥ 3 is that there exists a single large value s n ≥ s(1 -1 µ ), for which a large zone can be accommodated with low communication cost by the Square routine. When µ < 3 however, this guarantee is weakened to the existence of at most p = 1 µ-2 values whose sum is at least s(1 -1 µ ). It is thus necessary to design routines that can accommodate all p zones with low communication cost.

For NRRP, where we choose µ = 5 2 (and thus p = 2), we introduce the routine Tripartition, shown on Figure 8. This routine is used when (with the notations of Lemma 9) k = n -1, and divides the current rectangle R into two terminal zones Z 1 and Z 3 which are assigned to the two large remaining values, and one rectangle R 1 on which the recursive call is made.

Tripartition : The routine T ripartition(R, α, β) returns three rectangles R 1 , Z 2 and Z 3 of respective areas αs(R), βs(R) and (1-α-β)s(R) as depicted on Figure 8. In NRRP, the recursive call is thus only made on R 1 and we only need to prove ρ(R 1 ) ≤ 5 2 as an invariant on the aspect ratios of rectangles. 2) Approximation ratios: Changing the value of µ has also an effect on the approximation ratios for the individual zones returned by the algorithm. For the case where the algorithm returns a rectangle, the generalized version of Lemma 4 is stated in Lemma 10, with an identical proof.

Lemma 10. Let R be a rectangle. If ρ(R) ≤ µ then:

p(R) 2 s(R) ≤ µ + 1 2 √ µ .
This function is increasing for the interesting values of µ, and thus decreasing µ improves the approximation ratio in the case where all the produced zones are shaped as rectangles. However, as shown above, this is not the limiting worst-case; indeed, 2 √ 3 ( µ+1 2 √ µ for µ = 3) is smaller than 3 2 . Hence one could expect that slightly increasing this part of the bound with a larger value of µ would allow to lower the approximation ratios obtained for zones returned as complements of a square, given in Lemma 6. However the generalized version, Lemma 11, does not allow this.

Lemma 11. Let R be a rectangle such that ρ(R) ≤ µ, α ∈ [0, 1] and R , Z = Square(R, α). If α ≤ 1 µρ(R) then p(Z) 2 √ s(Z) ≤ max( µ µ-1 , √ µ(µ+1) 2 √ µ-1 ).
Proof. Let us denote ρ(Z) = ρ and s(Z) = s and let us suppose that h = h(R) ≤ w = w(R) without loss of generality. Note that the covering rectangle of Z is R and therefore w = w(Z) and h = h(Z).

In this case w = ρh, what implies p(Z) = (1 + ρ)h. By definition of Square, s = (1 -α)s(R) and s(R) = hw = ρh 2 . Therefore s = (1 -α)ρh 2 ≥ (1 -1 µρ(Z) )ρh 2 and hence p(Z) 2 S(Z) ≤ (1 + ρ)h 2 (1 -1/µρ)ρh 2 ≤ √ µ(1 + ρ) 2 (µρ -1) .
One can prove that the function

x -→ √ µ(1+x) 2 √ µx-1 is de- creasing on [1, 1 + 2 µ ] and increasing on [1 + 2 µ , µ] and then √ µ(1+x) 2 √ µx-1 ≤ max( µ µ-1 , √ µ(µ+1) 2 √ µ-1 ). Then p(Z) 2 √ S(Z) ≤ max( µ µ-1 , √ µ(µ+1) 2 √ µ-1 ). If µ < 3, then max( µ µ-1 , √ µ(µ+1) 2 √ µ-1 ) = µ µ-1 .
With such a value of µ, the dominant worst case is when the rectangle R is a perfect square (see Figure 9(a)) and the low value of µ implies the possibility for the square to take a large portion of the rectangle. If µ > 3, then max(

µ µ-1 , √ µ(µ+1) 2 √ µ-1 ) = √ µ(µ+1) 2 √ µ-1 .
With such a value of µ, the dominant worst case is when the rectangle R has an aspect ratio of µ (see Figure 9(b)) and the large value of µ implies a huge deformation of the rectangle what increases the approximation ratio.

1 √ µ 1 1 (a) 1 √ µ 1 µ (b)
Figure 9: Dominant worst case depending on µ for the Square operation.

In summary, x → x x-1 is decreasing on ]1, 3] and

x → √ x(x+1) 2 √
x-1 is increasing on [3, +∞[, therefore, with this proof technique and without a change in SN RRP , the value µ = 3 yields the best possible approximation ratio. It is thus necessary to include additional cases to improve over SN RRP .

3) Description of NRRP: Intuitively, choosing a larger value of µ appears to be counter productive (since it tends to produce taller and skinnier rectangles), and we thus work with µ = 5 2 , which is the lowest value of µ such that p = 2 in Lemma 9, allowing to use the Tripartition routine. However, as explained above, achieving an improved approximation ratio requires to design more sophisticated recursive partitionings to avoid the worst case which happens when the remaining large value is not large enough. The actual N RRP generates 10 sub-cases (see Figure 10), all of which are fully detailed in [START_REF] Beaumont | A new approximation algorithm for matrix partitioning in presence of strongly heterogeneous processors[END_REF]. The complete depiction of NRRP and of its subroutines is given in the supplementary material.

≥ 2 5 ≥ 2 5
Figure 10: All the cases in NRRP. In each case, the gray rectangle is the one on which the recursive call is made, the white ones are returned zones.

Theorem 12 (from [START_REF] Beaumont | A new approximation algorithm for matrix partitioning in presence of strongly heterogeneous processors[END_REF]). NRRP is a 2 √ 3 -approximation for PERI-SUM.

V. COMPARISON OF NRRP AND OPTIMAL

In this section, we compare the results returned by NRRP and the known optimal partitioning for a large number of computational power ratios (CP ratios), under two restrictions. First, we are concerned with the sum of half perimeters (SHP) only, which for convex rectilinear partitions is equivalent to the sum of the projections of partition areas along the axes, as described in Section II. Second, we compare NRRP to the optimal partitionings for two and three partition arrangements only. The first restriction is justified for several reasons. First, both NRRP and the known optimal partitionings deliver perfect load balance, leaving communication the only parameter to be optimized. Second, the SHP is the most simple yet useful metric for comparing the communication volume of two load-balanced partitionings. Finally, SHP can be used for any number of partitionings and thus is not limited to our second restriction. The second restriction is justified since it is only for two and three partition arrangements that the optimal partitionings are known for all CP ratios.

In order to carry out this comparison we implemented simulations that compute the SHP of partitionings returned by NRRP and the SHP of the corresponding optimal partitioning for a large number of CP ratios across all ratio ranges possible. This is done for the two and three partition cases separately. Due to the complexity of the highly-recursive NRRP algorithm, two simulation implementations were employed. The two implementations were written in Python and Java by separate team members who did not communicate until results were complete, after which results were compared, and any inconsistencies resolved.

A. Two Partitions

For two partitions there are only two optimal partitionings: the Straight Line partitioning, achieved with one call of the Guillotine subroutine of NRRP described in Section IV-B, and the Square Corner partitioning, achieved with one call of the Square subroutine of NRRP, also described in section IV-B. For CP ratios s Q : s P where s Q ≤ s P , normalized so that s Q + s P = 1, the Square Corner partitoning is optimal when s Q ≤ 0.25. The simulation analyzed CP ratios ranging from the only homogeneous case of s Q = s P = 0.5, to the most heterogeneous case of 0.005 : 0.995 in steps of 5.0 × 10 -3 . As NRRP can achieve both optimal partitionings in one step, calling subroutines that return known optimal shapes, it is not surprising that our simulations did not find any CP ratios where NRRP was not optimal.

B. Three Partitions

The three partition case is more complicated than the two partition case. There are three optimal partitionings as determined in Section III-B, but determining which one is optimal for a given ratio is not as straightforward as in the two processor case. With the same notations as in Section III, we denote the processors P, Q, R from the fastest to the slowest, and by s P , s Q , s R their respective speeds. Figure 11 shows the possible values of s R (on the x axis) and s Q (on the y axis) for CP ratios s R : s Q : s P where s R ≤ s Q ≤ s P , normalized so that s R + s Q + s P = 1. This space is dominated by a large region where the Block Rectangle partitioning is optimal, flanked by two regions along the y-axis (s R <≈ 0.086) where the Square Corner partitioning is optimal for values of s Q ≤ 0.25, and the Square Rectangle partitioning is optimal for values of s Q ≥ 0.25. At s Q = 0.25 both are optimal provided s R <≈ 0.018. These s R values are presented as approximate but the exact values are known; they are the solutions of the quadratic expressions that define the various regions.

For three partitions, our simulation analyzed 867 CP ratios ranging from the only homogeneous case of s Q = s R = s P = 0.333333 to the most heterogeneous case of 0.003333 : 0.003333 : 0.993333. We started with the homogeneous case and incremented s P by 0.01 until we reached s P = 0.993333. For each s P , we started with 1-s P 2 : 1-s P 2 : s P and then incremented s Q by 0.01 and decremented s R by 0.01, while s Q < s P (and s R > 0). This resulted in different numbers of CP ratios for each of the 67 values of s P , totaling different 867 CP ratios.

For each CP ratio, if the difference between the SHP returned by NRRP and the optimal SHP was greater than 2.0 × 10 -6 , then NRRP is considered to be non-optimal. We chose this value as all ratios were calculated to a precision of 1.0×10 -6 . Under these conditions NRRP is optimal for a total of 596 CP ratios (just over 2/3 of all CP ratios). Although nonoptimal for 276 CP ratios, NRRP is never far from optimal

< s R ≤ s Q ≤ s P and s R + s Q + s P = 1, s R ≤ 1 3 and s Q < 1 2 .
(average of 1.94%, maximum of 7.49% above optimal). This maximum is below the bound of The y-axis of Figure 12 shows the percent of all CP ratios where NRRP is not optimal (n = 867), and the x-axis shows how much above optimal these CP ratios are, from 0 to 7.5%, in intervals of 0.5%. The height of each bar indicates the total percent of all CP ratios where NRRP is not optimal for each interval. The sub-bars indicate the percent of CP ratios where each partitioning shape is optimal for each interval. For instance, the left-most bar in Figure 12 shows that just over 10% of all CP ratios where NRRP is not optimal are less than 0.5% above optimal. Of this, Square Rectangle is optimal for just under 8%, Square Corner is optimal for just over 1%, and Block Rectangle is optimal for less than 1%. We discuss the case of each shape in detail below.

• NRRP is not optimal for 69 CP ratios where Square Rectangle is the optimal shape. Square Rectangle is dominant where NRRP is very nearly optimal -within 0.5% (xaxis, Figure 12). Our simulations showed that these are CP ratios where s P is in the lower portion of its range, s Q is in the upper half of its range, and s R is extremely small (s P ≤ 0.733333, s Q ≥ 0.253333, s R ≤ 0.028333). • NRRP is not optimal for 112 CP ratios where Square Corner is the optimal shape. Square Corner is dominant when NRRP is between 0.5% and 3.0% from optimal (xaxis, Figure 12). Our simulations showed that these are CP ratios where s P is above its Square Rectangle range, s Q is in the lower half of its range, and s R is small (s P ≥ 0.743333, s Q ≤ 0.248333, s R ≤ 0.083333). • NRRP is not optimal for 95 CP ratios where Block Rectangle is the optimal shape. Block Rectangle is dominant when NRRP is between 3.0% and 7.5% from optimal (x-axis, Figure 12). Our simulations showed that these are CP ratios where s P is in the middle of its range, s Q varies considerably, and s R is small (0.453333 ≤ s P ≤ 0.823333, 0.088333 ≤ s Q ≤ 0.458333, s R ≤ 0.123333). Although there is some overlap in CP ratio ranges where the three optimal shapes are dominant, some conclusions can be drawn. For instance, NRRP is always optimal when s R > 0.123333. Additionally, the values of s Q and s P determine when Square Corner or Square Rectangle are optimal. Square Rectangle is only a possible optimal shape when s P <≈ 0.74 and s Q >≈ 0.25. Square-Corner is only a possible optimal shape when s P >≈ 0.74 and s Q <≈ 0.25. This 'y-boundary' is also clearly seen in Figure 11, which shows the optimal shape for all possible ratios, regardless of whether NRRP is optimal or not.

C. More Than 3 Partitions

We do not investigate cases of 4 (or more) partitions as it is not known how many optimal partitioning shapes there are, what they are, and for what ratios they are optimal. In our experience the three partition case is much more complex than the two partition case. We have no reason to believe that the four partition case will be any different, and thus it is beyond the scope of this work. Additionally, it is difficult to speculate as to how NRRP will perform in this case, but it is possible that NRRP performance decreases with increasing partitioning complexity as for two partitions we saw that NRRP is optimal for all ratios, but for three partitions we saw that NRRP was not optimal for approximately 1 3 of ratios. Nonetheless, all NRRP partitionings will be within 2 √ 3 of the optimal. This is one of the most important strengths of the NRRP algorithm -even where the optimal partitionings are not known (which currently is all partitionings with more than three partitions), NRRP returns a partitioning which is no more than a factor of 2 √ 3 of the optimal. However, in [START_REF] Beaumont | A new approximation algorithm for matrix partitioning in presence of strongly heterogeneous processors[END_REF], some simulations with numbers of processors greater than 3 has been performed (up to 64 processors). In this set of experiments, three kinds of processors have been considered, each with different speed (one representing CPU, another GPU and the last accelerators) and for each value of n (number of processors), different repartitions of CPUs/GPUs/accelerators have been tested. As stated previously, for such numbers of processors, optimal solutions are not known and the comparison relies on the lower bound in Equation 1. With such settings, NRRP returns, on average, solutions whose SHP is within 5% of the lower bound. The maximum ratio between a solution returned by NRRP and the bound is 1.106. Note that solutions of RRP (Algorithm 1) have also been considered, with behavior similar to the one of NRRP, except that the worst case scenario is significantly worse (around 1.3).

VI. CONCLUSION

The problem of partitioning a matrix into a set of submatrices to address simultaneously the problems of load balancing and communication minimization for data parallel applications on heterogeneous platforms has been proved to be NP-Complete. In this paper, we present recent approaches Figure 12: Percent from optimal for all computational power ratios where NRRP is not optimal, and which partition shape is optimal in these cases.

that relax the restriction that all partitions should be shaped as rectangles. The first approach uses an original mathematical technique to find the exact optimal partitioning for a small number of processors. This is important as in our experience there is an appetite for exact solutions, even in this restricted setting, particularly as it is applicable to current CPU/GPU/accelerator configurations. Due to the complexity of this technique, optimal algorithms are known for up to three heterogeneous resources only. However, even at this small scale, we can observe that optimal partitions are often non-rectangular and sometimes non-intuitive. This justifies the extension of the non-rectangular partitioning problem to find approximate solutions for any number of resources.

In this direction, we propose two algorithms, i.e. Simple Non-Rectangular Recursive Partitioning (SNRPP), which is a 3 2 approximation algorithm, and Non-Rectangular Recursive Partitioning (NRPP), which is a 2 √ 3 approximation algorithm. While sub-optimal, these algorithms return results for arbitrary numbers of partitions, an advantage for problems with larger number of processors where the shape of optimal solutions is unknown. We use the first exact approach to analyze how close to the known optimal solutions the NRRP algorithm is for small numbers of partitions. We show that this approximation algorithm produces partitions that are often very close to optimal. Our improved algorithms have the potential for significant practical impact: in many applications, matrix multiplication operations are used iteratively a large number of times, so that even small improvements have a significant impact. We have observed this impact in our work, where a small difference in theory leads to a large difference in practice, particularly when problem size and network topology are taken into account. In one example, at a computational power ratio where the optimal partition shape outperforms the NRRP partition by only 1.86%, the optimal shape outperforms the NRRP partition by 23.63%. Thus, a nearly optimal solution in terms of communication cost can turn out to be far from the optimal in terms of real-world execution time.

This work opens several directions of future research. As far as optimal partitionings are concerned, it would be interesting to study cases with a larger number of resources (4, 5, ...) and to study the actual complexity of finding an optimal partition for the general case. Improvements of the approximation ratio might also be possible, for example by considering smaller values for the parameter µ. Another direction for future research would be to study other related models, in particular to consider the case of sparse rather than dense matrices, or to consider more precise communication models. In the sparse matrices case, the processing cost for a processor is related not to the area assigned to it, but to the number of non-zero elements in its zone, which introduces a new source of heterogeneity related to input matrices. More precise communication models would require network topology to be taken into account explicitly, as proposed in [START_REF] Deflumere | Optimal partitioning for parallel matrix computation on a small number of abstract heterogeneous processors[END_REF]. In this case, similar to those studied here, non-intuitive shapes are more dominant as heterogeneity increases, which provides a strong motivation to continue working on finding optimal solutions, even in the case of a small number of heterogeneous resources, but with explicit topologies.

Finally, the natural next step is to provide practical implementations of these algorithms in actual matrix multiplication kernels. This work has already been started [START_REF] Lambert | On the Effect of Replication of Input Files on the Efficiency and the Robustness of a Set of Computations[END_REF], with work in two directions. First, this requires a move from continuous partitionings to the discrete partitioning of a [1, N ] × [1, N ] matrix. To the best of our knowledge, finding approximation algorithms for the discrete problem has not been addressed yet. Second, practical implementations also require dynamic work-stealing strategies to correct for the unpredictability of computing times inherent to large-scale machines.

Figure 2 :

 2 Figure 2: The candidate partition shapes identified as potentially optimal three-processor shapes in [12]. (1) Square Corner (2) Rectangle Corner (3) Square Rectangle (4) Block Rectangle (5) L Rectangle (6) Straight Line

Figure 3 :Figure 4 :

 34 Figure3: The optimal data partitions for a rectangle with three processors.

  then the Square Corner solution of Figure 3 is feasible, and has cost exactly a + b + 2 √ s Q + 2 √ s R . Since Q and R are the two smallest size values, this solution dominates S. (iv) Exactly 2 resources T and U hold a complete side in S.

Figure 3

 3 Figure 3 is feasible, and has cost 2a + b + 2 √ s R . Since s V ≥ s S , S is dominated. (v) All 3 resources hold a complete side in S. If they hold a complete column, then C ≥ 3a + b, otherwise C ≥ 3b + a ≥ 3a + b. Then, S is dominated by the Straight Line solution.

Figure 5 :

 5 Figure 5: Two partitions of a square with the instance { , 1-}.

Figure 6 :

 6 Figure 6: An illustration of the Guillotine (left) and Square (right) routines.

Figure 7 :

 7 Figure 7: Worst case scenario for SNRRP.

R 1 Z 2 Z 3 RFigure 8 :

 1238 Figure 8: An illustration of the Tripartition routine.

Figure 11 :

 11 Figure 11: Regions where Square Corner, Square Rectangle, and Block Rectangle partitionings are optimal for three partitions, with x = s R and y = s Q . Note that as 0< s R ≤ s Q ≤ s P and s R + s Q + s P = 1, s R ≤ 1 3 and s Q < 1 2 .

2

 2

√ 3

 3 established by Theorem 18 by a factor of over 2.

VII. ACKNOWLEDGEMENTS This publication has emanated from research conducted with the financial support of Science Foundation Ireland (SFI) under Grant Number 14/IA/2474. This work is partially supported by EU under the COST Program Action IC1305: Network for Sustainable Ultrascale Computing (NESUS).