Alexey Brodoline
email: alexey.brodoline@umontpellier.fr

Nitin Rawat

Daniel Alexandre

Nicolas Cubedo

Michel Gross

4D compressive sensing holographic imaging of small moving objects with multiple illuminations

4D compressive sensing holographic imaging of small moving objects with multiple illuminations

Introduction

Digital holography (DH) enables to retrieve the amplitude and the phase of the optical field scattered by an object from the interference between this field and a reference beam [START_REF] Gabor | A new microscopic principle[END_REF]. The resulting interferogram, called the hologram, is recorded on an image sensor, making possible numerical reconstruction of the complex field [START_REF] Schnars | Direct recording of holograms by a ccd target and numerical reconstruction[END_REF].

For this purpose, the field on the sensor is backpropagated within a series of planes adjacent to the object. In this way, a spatial distribution of the field scattered in the vicinity of the object is obtained. Nevertheless, this 3D field does not fit the 3D distribution of the scattering density of the object.

As was shown in [START_REF] Atlan | Heterodyne holographic microscopy of gold particles[END_REF], a gold nanoparticle imaged by digital holographic microscopy (DHM) will give a cone of light after field backpropagation. The angle of the cone corresponds to the numerical aperture of the microscope objective, while its apex is located at the position of the gold bead. Holography, which is able to localize this apex, is thus well suited to track small particles, and was used for this purpose in a large number of works, especially flow imaging [START_REF] Buraga-Lefebvre | Application of wavelet transform to hologram analysis: three-dimensional location of particles[END_REF][START_REF] Soulez | Inverse-problem approach for particle digital holography: accurate location based on local optimization[END_REF][START_REF] Tian | Quantitative measurement of size and three-dimensional position of fastmoving bubbles in air-water mixture flows using digital holography[END_REF][START_REF] Memmolo | Twin-beams digital holography for 3d tracking and quantitative phasecontrast microscopy in microfluidics[END_REF].

Due to its phase retrieval capability, holography is a valuable complement to standard microscopic techniques and was widely applied to imaging of biological samples [START_REF] Cuche | Simultaneous amplitude-contrast and quantitative phase-contrast microscopy by numerical reconstruction of fresnel off-axis holograms[END_REF][START_REF] Kang | Combined confocal raman and quantitative phase microscopy system for biomedical diagnosis[END_REF][START_REF] Quan | Multimodal digital holographic microscopy for widefield fluorescence and 3d phase imaging[END_REF]. Quantitative phase imaging techniques are able to provide an estimate of the thickness or of the dry mass of a living cell [START_REF] Rappaz | Simultaneous cell morphometry and refractive index measurement with dualwavelength digital holographic microscopy and dye-enhanced dispersion of perfusion medium[END_REF][START_REF] Kemper | Digital holographic microscopy for live cell applications and technical inspection[END_REF][START_REF] Popescu | Optical imaging of cell mass and growth dynamics[END_REF].

If a true 3D reconstruction of the object is required, techniques such as diffraction tomography should be considered. By recording a series of holograms for different directions of illumination and by performing a tomographic reconstruction, the 3D optical index distribution of the object can be computed [START_REF] Debailleul | High-resolution threedimensional tomographic diffractive microscopy of transparent inorganic and biological samples[END_REF][START_REF] Sung | Optical diffraction tomography for high resolution live cell imaging[END_REF]. However, because it requires to record a large number of holograms (≥ 100), tomography is not well suited for the imaging of moving objects.

In the case where there is a basis in which the 3D scattering density is sparse enough, it is also possible to reconstruct a 3D image of the object from a single hologram by compressive sensing (CS) [START_REF] Donoho | Compressed sensing[END_REF][START_REF] Candes | Stable signal recovery from incomplete and inaccurate measurements[END_REF]. This point has been demonstrated by Brady et al. [START_REF] Brady | Compressive holography[END_REF] who have imaged by holography an immobile object whose representation is sparse in the wavelet basis. Since this work, DH coupled to the CS has given rise to many developments. We can mention compressive in-line holography [START_REF] Denis | Inline hologram reconstruction with sparsity constraints[END_REF], compressive Fresnel holography [START_REF] Rivenson | Compressive fresnel holography[END_REF], off-axis frequency shifting holography [START_REF] Marim | Compressed sensing with offaxis frequency-shifting holography[END_REF], off-axis holography of diffuse objects [START_REF] Choi | Compressive holography of diffuse objects[END_REF], video-rate microscopic tomography [START_REF] Hahn | Video-rate compressive holographic microscopic tomography[END_REF], or incoherent holographic microscopy [START_REF] Cossairt | Compressive reconstruction for 3d incoherent holographic microscopy[END_REF]. Recently, a combination of compressive holography and coded exposure enabled the reconstruction of 3D videos of small sparse objects [START_REF] Wang | Compressive holographic video[END_REF].

In this work, we focused on 3D imaging of blood circulation in zebrafish (Danio rerio) larvae, which implies to localize in space and in time the positions of the red blood cells (RBCs). Due to the fact that their refractive index is close to the one of the surrounding medium, the cone of the scattered light is narrow. For this reason, the RBCs can not be precisely localized with classical holographic backpropagation. In [START_REF] Donnarumma | 4d holographic microscopy of zebrafish larvae microcirculation[END_REF], we proposed a technique based on a specific holographic setup and a "cleaning" algorithm, which enabled to compute the positions of the RBCs in 3D, as well as the shape of the blood vessels. Nevertheless, the drawback of this algorithm was the relatively long computing time. In [START_REF] Brodoline | 4d compressive sensing holographic microscopy imaging of small moving objects[END_REF], we described an algorithm which enabled to speed up calculations by a factor of ∼ 500, based on compressive sensing.

Indeed, RBCs can be viewed as small scattering objects that move within a larger motionless object (the anesthetized zebrafish larva). In the basis of the direct 3D space (i.e. x, y, z), their 3D scattering density is sparse, which allows to apply compressive sensing reconstruction methods. In this specific case, the use of sparsifying operators [START_REF] Bettens | Studies on the sparsifying operator in compressive digital holography[END_REF] is not required. We used an algorithm based on orthogonal matching pursuit (OMP). The case of objects composed of point sources was also studied by Liu et al. [START_REF] Liu | Compressive holography algorithm for the objects composed of point sources[END_REF], where a more complex reconstruction algorithm was proposed based on similar compact sensing matrix computation.

By performing calculations from a small number of holograms, we were able to obtain the 3D positions of the RBCs at one camera frame. By performing a series of reconstructions from a long sequence of holograms, we got movies of the 3D positions of RBCs, and by averaging in time, we determined the shape of the blood vessels from the trajectories of the RBCs.

In this paper, we will give further details on the experimental configuration and the reconstruction procedure. We will, in particular, describe a method for the localization of the RBCs that uses the fact that the object is illuminated by several beams, which results in a more robust reconstruction algorithm.

Experimental setup

The experiment is carried out in an off-axis configuration similar to previous work [START_REF] Donnarumma | 4d holographic microscopy of zebrafish larvae microcirculation[END_REF]. To increase the angular diversity of illumination, the sample was imaged with a higher NA microscope objective (20x/NA=0.5 water immersion, Zeiss) and with three illumination beams, instead of NA=0.3 and two illuminations in [START_REF] Donnarumma | 4d holographic microscopy of zebrafish larvae microcirculation[END_REF]. The holographic signal was acquired by a Mikrotron EoSens CL camera (1280 × 1024 pixels, 14 µm, 12 Bits, 200 frames/s). The laser diode power at the sample was less than 1 mW at λ = 662 nm.

The experimental setup is depicted in Fig. 1. The illumination beam E IL generated by a laser diode is divided into three using beam splitters (BS) and mirrors (M). The object is thus illuminated from three different directions. The intensities of the beams can be equalized with neutral density filters. The scattered light E is collected by the microscope objective (MO) and imaged on the camera, in plane C , where it interferes with the reference field E R . The combining BS is tilted with an angle θ to perform off-axis holography. The camera records the resulting interferogram

I = |E C + E R | 2
, where E C is the field scattered by the object in the camera plane C .

We have validated the technique by imaging the microcirculation in zebrafish (Danio rerio) larvae. Unpigmented larvae were obtained from spontaneous spawning of casper young adult fish pairs and grown at 28.5 • C. They were anesthetized with tricaine at five days and mounted in lateral view in 1% low melting point agarose in a 35 mm glass ibidi bottom petri dish. The trunk blood flow was imaged immediately. Animals were handled according to standard procedures [START_REF] Westerfield | The zebrafish book : a guide for the laboratory use of zebrafish (Danio rerio)[END_REF].

Selection of the moving objects

The zebrafish larva is an optically thick object composed of tissues, organs and bones scattering light, which results in a signal irrelevant for our applications. In order to select the moving RBCs and eliminate this background, we considered sequences of interferograms I m with m = 0, 1, 2... corresponding to successive camera frames. The resulting hologram H C in the camera plane C is then obtained by combining the I m with coefficients whose sum is zero:

H C (x, y) = m η m I m (x, y) with m η m = 0. (1)
For example,

H C = (I m + I m-1 ) -(I m-2 + I m-3 ), and
η m = (1, 1, -1, -1).
In this way, the immobile background is subtracted, and the RBCs whose positions changed across the frames remain in the hologram. The number of I m and the interval between them can be changed in order to increase the signal amplitude. This approach is very similar to the one of Schwerte and Pelster [START_REF] Schwerte | Digital motion analysis as a tool for analysing the shape and performance of the circulatory system in transparent animals[END_REF] who highlighted blood circulation in zebrafish larvae by calculating differences of successive frames in white light video movies. This approach was also used with holography in [START_REF] Gao | In vivo imaging and quantitative analysis of zebrafish embryos by digital holographic microscopy[END_REF].

In practice, we considered combinations of six holograms, according to the following equation (with N = 6):

H C (x, y) = N -1 k=0 sin( 2πk N ) I m (x, y). (2)
If we consider a higher N , we will average on more frames, which will result in more signal but can as well increase the noise.

Formulation of the compressive sensing equations

A simplified outline of the setup with a single illumination is depicted in Fig. 2. The imaged objects are located in a volume V and illuminated by a laser diode L. The scattered light E is collected by the MO and imaged on the camera, in plane C . It interferes with the reference field E R to form the hologram

I = |E C + E R | 2
, where E C is the field scattered by the objects in the camera plane. Since the experiment is carried out in off-axis geometry, by filtering in the Fourier space [START_REF] Cuche | Spatial filtering for zero-order and twin-image elimination in digital off-axis holography[END_REF], we can select the +1 order component of the hologram, corresponding to the

E C E * R term.
By considering that the reference field is a flat field (|E R | 1), we thus calculate the field E C and eliminate the phase effects related to eventual aberrations of the MO [START_REF] Colomb | Automatic procedure for aberration compensation in digital holographic microscopy and applications to specimen shape compensation[END_REF], to the off-axis configuration and to the curvature of the reference wave [START_REF] Verrier | Holographic microscopy reconstruction in both object and image half-spaces with an undistorted three-dimensional grid[END_REF]. We get therefore the field E C (x, y) in the plane C, conjugate of the plane C by the MO. We get as well the field ẼC in the Fourier space, ẼC (k

x , k y ) = FT [ E C (x, y) ],
where FT is the discrete Fourier transform operator.

From now, we can combine several holograms in order to eliminate the background and select only the moving RBCs. In the following, the fields E C , E C and their Fourier transforms designate the fields scattered by the moving objects only.

The RBCs that scatter light, act like small sources of field S(x, y, z) located in volume V [Fig. 2]. According to angular spectrum propagation equations [START_REF] Yu | Wavelength-scanning digital interference holography for tomographic three-dimensional imaging by use of the angular spectrum method[END_REF], the field E C (x, y) radiated by the sources S in the plane C at z=0 is E C (x, y) = z FT -1 exp(+ik z .z) FT [ S(x, y, z) ], (where FT -1 is the inverse Fourier transform operator). With the linearity of the Fourier transform, we get the expression of the scattered field spectrum in the camera conjugate plane:

ẼC (k x , k y ) = z exp(+ik z .z) FT [ S(x, y, z) ] , (3)
where k z is the z component of the wavevector k, i.e.,

k z = k 2 -k 2 x -k 2 y with k = |k| = 2πn m /λ
, where λ is the wavelength in vacuum and n m 1.33 the refractive index of the medium.

Equation (3) can be written as a linear system of equations that connects the 3D field scattered by the objects to the field spectrum in the camera conjugate plane:

u = Av, (4)
where v ≡ S(x, y, z) and u ≡ ẼC (k x , k y ) are vectors of dimensions 768 × 768 × 160 and 768 × 768, respectively, and A is a matrix of dimension 768 4 × 160 [see Fig. 3(a)]. This system is strongly under-determined since the dimensions of vectors u and v correspond to the number of pixels of the calculation grid in plane C and the number of voxels of the volume V, respectively. The compression ratio (CR) between the vectors is large, since it corresponds to the number of planes z considered for the 3D reconstruction (i.e., CR= 160).

We should note that the matrix A is too large to be calculated. To compute u = Av, (3) must be used. Similarly, to compute v = A † u, we must use (5) for each z:

S(x, y, z) = FT -1 exp(-ik z .z) ẼC (k x , k y ) . (5)
The moving RBCs represent a very small part of the imaged volume, and the vector v is highly sparse in the basis of the positions x, y, z. It is then possible to solve (4) iteratively using compressive sensing.

Orthogonal matching pursuit algorithm

The CS method we used here is based on the orthogonal matching pursuit (OMP) [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF]. It is described in detail in our previous article [START_REF] Brodoline | 4d compressive sensing holographic microscopy imaging of small moving objects[END_REF]. Let us recall the main steps of the algorithm. It can be schematized by the equation:

u (n+1) = u (n) -a n AT (n) A † u (n) , (6)
where n is the iteration index. The initialization u (1) = u ≡ ẼC is done with the object field in the Fourier space obtained from the +1 order of the recorded hologram after phase corrections and background subtraction.

The algorithm calculates iteratively the 3D sources v (n) = A † u (n) , and the 2D field Av (n) radiated by the sources, by selecting within the sources v (n) the brightest ones with the operator T (n) . The 2D field radiated by the selected sources (i.e. AT (n) A † u (n) ) is subtracted from the 2D field u (n) in order to get the 2D field for the next iteration u (n+1) . In this way, the 2D hologram in Fourier space is cleaned at each iteration, while the positions of the selected sources are saved.

We should remark that in order to optimize and stabilize the cleaning procedure, the selected sources T (n) v (n) are multiplied by a normalization factor a n . It is adjusted so that the resulting vector u (n+1) remains orthogonal to the vector u n = AT (n) A † u (n) radiated by the selected sources [Fig. 3(b)]. The energy of the 2D field |u (n+1) | 2 is thus minimized after cleaning.

To achieve the condition of orthogonality, we have to satisfy (u (n) -a n u n ) u n = 0, where is the scalar product. We get

a n = u n u (n) |u n | 2 . ( 7
)
At each iteration step n, the total energy (summed over all pixels) remaining in the 2D field |u (n) | 2 decreases. When this energy becomes sufficiently low, e.g., when it represents only 20% of the initial energy, the iterative process stops. The final 3D distribution v of the sources is given by the summation of the sources selected at each iteration:

v = n a n T (n) A † u (n) . (8)
It represents the 3D distribution of the RBCs at one camera frame.

6 Operator T (n) for the selection of the sources

The reconstruction speed and quality depend of the operator T (n) that selects at each iteration n the brightest sources. There are several ways to define this operator. In [START_REF] Brodoline | 4d compressive sensing holographic microscopy imaging of small moving objects[END_REF], the criterion to select the RBCs was the maximum of the energy of the complex field.

A second manner to proceed is to use the energy correlation between the three laser beams illuminating the object. Let us give here more details on these methods of selection of the RBCs.

Selection of the sources based on energy

In order to find the positions of the RBCs, we can calculate the energy of the complex 3D field and look for its maxima. Indeed, the RBCs are located at the intersections of apexes of the cones of scattered light there energy is maximum. We should remark that the localization is possible only due to the multiple illumination. Since the refractive index of RBCs is very close to the one of the surrounding medium, the angle of the scattering cones is small, and thus the localization with a single beam would not be precise enough.

Let us consider the Fourier transform of the off-axis hologram, represented on a 1280 × 1280 calculation grid [Fig. 4(a)]. The zero grating order is on the center, while the +1 and -1 orders that represent the rear pupil of the microscope objective are on the top and on the bottom, respectively. The three bright spots inside the pupil correspond to the ballistic light of the illumination beams. Note that the objective has a phase ring which can be seen on the figure.

Since it does not overlap with the beams, it has no effect on our reconstructions.

The +1 grating order signal (which corresponds to ẼC (k x , k y )) was selected by cutting a disk of radius 340 pixels within the Fourier space. This disk was copied on a 768 × 768 calculation grid [Fig. 4 (b)]. The zero order was filtered by a circular mask of radius 170 to avoid signal mixing. The center of the pupil with the phase ring was as well removed with a disk of radius 195 pixels.

The algorithm starts with the filtered +1 order u (1) ≡ ẼC (k x , k y ). At each iteration step n, after the calculation of the 3D field v (n) = A † u (n) with ( 5), the energy of the sources |v

(n) | 2 ≡ |S (n) (x, y, z)| 2 is calculated for each voxel (x, y, z). Then a threshold of energy E (n) th is chosen, e.g. E (n) th = 0.8 E (n) max , where E (n) max is the maximum of the vector |v (n) | 2 .
Afterwards, the operator

T (n) E
selects the voxels (x , y , z ) whose energy is above the threshold E (n) th and adds to these voxels the voxels located within a predefined radius (e.g., r RBC = 4 pixels) in the same z = z plane. The selected sources T

(n) E A † u (n) ≡ S n (x, y, z) are thus • S n (x, y, z) = S (n) (x, y, z) if |S (n) (x , y , z )| 2 > E (n)
th , (x -x) 2 + (y -y) 2 ≤ r RBC , and z=z ,

• S n (x, y, z) = 0 otherwise.

Two criteria can be used to terminate the algorithm. In [START_REF] Brodoline | 4d compressive sensing holographic microscopy imaging of small moving objects[END_REF], the stopping criterion was the value of the total energy of the hologram. When it becomes lower than 20% of the initial energy, i.e., if |u (n) | 2 < 0.2 |u (1) | 2 , the algorithm ends. Here, we considered the ratio between the maximum energy at the current iteration and the maximum at the first iteration. When it becomes sufficiently low, e.g., if

E (n) max /E (1)
max < 0.06, the calculation will stop. In practice, the second criterion is reached before the first. The final 3D sources v ≡ S(x, y, z) are then obtained by a summation of all the detected sources and by taking into account the normalization coefficients:

S(x, y, z) = n a n S n (x, y, z). (9)

Selection of the sources with multiple illuminations

A second approach is to consider separately the three laser beams illuminating the object from different directions. The separation can be performed in the Fourier space, while the +1 order of the hologram is selected. Figure 4(c) shows the crop of the +1 order, where three circular masks (in red, green and blue) are applied to select the light scattered by the illumination beams. The bright peaks correspond to ballistic light, while the surrounding halo constitutes scattered light. The complex 3D field can then be computed for the three beams separately, giving the vectors v (n,i) ≡ S (n,i) (x, y, z) with i = 1, 2, 3. The intensities of the three fields in the direct space in plane C are represented in Fig. 4(d).

The principle of the OMP algorithm is the same as previously. Nevertheless, the transforms associated with A and A † are applied to the three fields separately, and the reconstruction of the vectors u (n,i) (with i = 1, 2, 3) in the Fourier space with (3) has to take into account the selection masks.

Another major change is the criterion of selection of the sources S n (x, y, z) with the operator T (n) C . After the calculation of the vectors S (n,i) in the object space, we calculate for each voxel the product of the energies of the three fields, which we call the correlation Corr 3 , similar to what was done in [START_REF] Donnarumma | 4d holographic microscopy of zebrafish larvae microcirculation[END_REF]:

Corr 3 n (x, y, z) = |S (n,1) (x, y, z)| 2 (10) × |S (n,2) (x, y, z)| 2 × |S (n,3) (x, y, z)| 2 .
Since we have a product of three energies, to get an equivalency with the previous method, we should consider the cube root Corr n = 3 Corr 3 n at each pixel. We define then a threshold C th and adds to these voxels the voxels located within a predefined radius (e.g., r RBC = 4 pixels) in the same z = z plane. The selected sources T

(n) C A † u (n) ≡ S n (x, y, z) are thus • S n (x, y, z)= S (n,i) (x, y, z) i if Corr n (x , y , z ) > C (n)
th , (x -x) 2 + (y -y) 2 ≤ r RBC , and z=z ,

• S n (x, y, z) = 0 otherwise, where . i represents the average among the three fields.

The remaining calculations such as the cleaning of u (n) with the use of the coefficient a n and the summation of all the sources after the last iteration are the same as for the energy. The CS calculations were conducted on a 768 × 768 × 160 tridimensional grid and on 256 camera frames in approximately 2 h using GPU (NVidia GTX Titan Xp) and CUDA libraries (cuFFT, cuBlas...). Voxel size was 0.6 µm in x and y, and 1.8 µm in z. The 3D visualizations were performed with a homemade CUDA program based on an NVidia code example (VolumeRender). Figure 5(a) shows the 3D positions of the RBCs at one camera frame obtained from a combination of six holograms modulated by a sinus and using the operator

T (n) C
that selects the brightest sources taking account of the correlation among the three illumination beams. It can be noticed that each RBC produces two bright spots due to the image subtraction. A complete 360 • rotation can be seen in (Visualization 1).

Figures 5(b) and 5(c) show the 3D images of the shape of the blood vessels, obtained by averaging the positions of the RBCs over 256 camera frames. The reconstructed vascular system compares favorably to images obtained by confocal microscopy [START_REF] Isogai | The vascular anatomy of the developing zebrafish: an atlas of embryonic and early larval development[END_REF]. We should remark that we do not reconstruct the blood vessels walls like it can be done in fluorescence microscopy. Nonetheless, the shape of the main trunk vessels is well reconstructed. A complete 360 • rotation of the vascular system can be found in (Visualization 2).

In order to visualize the 3D motion of the RBCs inside the blood vessels, we superimposed the nonaveraged and averaged reconstructions [Fig. 5(d)]. A 360 • rotation of this sequence of 3D images is shown in Visualization 3. In the movie, the direction of blood flow in the dorsal aorta, caudal vein, intersegmental vessels and dorsal longitudinal anastomotic vessels can be determined. Since the camera has a relatively high frame rate (200 images/s), the RBCs motion appears fluid. reconstruction that considers the three illumination beams [Fig. 6(a)]. We can remark that the difference is very small. Nevertheless, the image obtained with T (n) E seems to be noisier, with artifacts along the field propagation direction z (white arrows). considered. Nevertheless, the computing time of one iteration is about twice longer, and therefore the total calculation times remain equivalent for both cases. (n) th = 0.5 Corr n,max ). In this case, the algorithm converges very fast, since an important amount of energy is removed at each iteration. Nevertheless, this reduces the quality of the 3D reconstructions. In Fig. 7(b), in the case where the reconstruction is based only on energy, the ratio variation is not stable and increases at some iterations. The second algorithm is more robust and does not present such non-uniformity.

Convergence of the algorithm

Figure 8 shows the reconstructions obtained with low thresholds (i.e., 0.5). In the case where only the energy was considered, the 3D image is incomplete. The dorsal vessels where the blood flow was not significant are not reconstructed. High-energy artifacts elongated in the z direction are as well present (white arrows). The reconstruction with the three illumination beams [Fig. 8(b)] is more robust, and the shape of the vascular system is nonetheless retrieved. In some frames however, the least bright RBCs in the dorsal vessels are not reconstructed, which results in a blinking image along the frames.

Conclusion

In this paper, we presented a reconstruction procedure based on CS and a specific DHM setup using a triple illumination of the object. This technique is able to efficiently retrieve the 3D positions of the RBCs in the zebrafish larva vascular system, since they represent small scattering objects moving in a large motionless volume. The reconstruction of the 3D positions of the RBCs as well as the shape of the perfused vascular system enables to perform in-vivo imaging of microcirculation.

We compared two methods of selection of the RBCs; one is based on the maximum of the energy of the hologram, while the second considers the correlations in energy between the three illumination beams. The last one takes explicitly into account the fact that the object is illuminated by three beams. It appears to be more robust to ill choices of threshold parame-ters and gives less noisy results.

We should remark that several parameters influence the amplitude of the reconstruction noise. For example, in order to eliminate the background, we considered combinations of six holograms with sinus coefficients η m . By combining more images, we will average on more frames and get therefore more signal with, at the same time, the risk of increasing noise. We could as well consider complex exponential coefficients η m . In this case, the RBCs would appear more elongated. The influence of the stopping criterion should also be taken into account. If the chosen threshold is too small, irrelevant low energy signal will be selected, resulting in reconstruction artifacts in the image.

The basis used for CS calculations is the direct 3D space (x, y, z) and since the object is sparse in this basis, the OMP method works quite well. The dictionary is simply composed of point sources at every (x, y, z) coordinate. Since the object is illuminated by several beams, the reconstruction does not require an orthogonal basis where every pair of voxels would have uncorrelated holograms. An optimized basis could be learned, but it would necessitate a large number of test holograms to constitute the dictionary as in [START_REF] Liu | Efficient spacetime sampling with pixel-wise coded exposure for high-speed imaging[END_REF][START_REF] Wang | Katsaggelos, and Oliver Cossairt. Dictionary-based phase retrieval for space-time super resolution using lens-free on-chip holographic video[END_REF].

We could also try more complex reconstruction algorithms on our data, in particular the one proposed by Liu et al. [START_REF] Liu | Compressive holography algorithm for the objects composed of point sources[END_REF], or regularized algorithms such as LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], thresholding algorithms such as FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], and bayesian algorithms [START_REF] Ji | Bayesian compressive sensing[END_REF]. We have nevertheless to take into account memory issues, since the dimensions of our sensing matrix are very large.

Figure 1 :

 1 Figure 1: Triple illumination off-axis holographic microscopy setup. BS, beam splitters; M, mirrors; MO, microscope objective; θ, off-axis angle; C , camera plane; E IL , illumination beam; E, scattered object field; and E R , reference field.

Figure 2 :

 2 Figure 2: Typical off-axis holographic microscopy setup. L, laser; C', camera; BS, beam splitter; M, mirror; MO, microscope objective; C, conjugate plane of plane C' by objective MO; V, imaged volume; S, a scattering object; E, object field; and E R , reference field.

Figure 3 :

 3 Figure 3: (a) Schematic representation of the transformations. Matrix A describes the emission by the 3D sources v ≡ S resulting in a 2D field in Fourier space u ≡ ẼC . Matrix A † represents the 3D field reconstruction from the 2D field in the camera conjugate plane. (b) Normalization procedure diagram: u (n) is the hologram at the iteration n, u n the field radiated by the sources selected at iteration n, a n the normalization coefficient, and u (n+1) the hologram at the iteration n + 1, orthogonal to u n .

Figure 4 :

 4 Figure 4: (a) Fourier transform of the off-axis hologram. The bright spots inside the +1 and -1 orders correspond to the three illumination beams. (b) Crop of the +1 order of the hologram, filtering of the center of the pupil and of the zero order. (c) Separation of the three illumination beams inside the +1 order with circular masks (red, green and blue). (d) Intensity of the three fields (red, green and blue) reconstructed in the camera conjugate plane.

  on the maximum value of Corr n , e.g., C (n) th = 0.8 Corr n,max . The operator T (n) C selects the voxels (x , y , z ) whose correlation is above the threshold C (n)

  images of zebrafish microcirculation

Figure 6 (

 6 b) shows the reconstruction obtained considering only the energy of the hologram (i.e., with the operator T (n) E ). It can be compared with the

Figure 5 :

 5 Figure 5: (a) Representation of RBCs considering a single 3D reconstruction at one camera frame. View at 0 • . (b,c) Representations of the vascular system considering a summation of a sequence of 256 holograms. Views at 96 • and 120 • . (d) The RBCs positions are superimposed with the shape of the blood vessels. Dorsal aorta (DA), caudal vein (CV), intersegmental vessels (ISV) and dorsal longitudinal anastomotic vessels (DLAV) can be clearly seen. View at 0 • . Reconstructions performed with operator T (n) C

Figure 7 (C

 7 Figure 7(a) shows the mean decrease in total energy of the hologram with time for operators T (n) E and T (n)C . The algorithm converges in fewer iterations when the correlation among the three illuminations is

Figure 6 :

 6 Figure 6: Comparison between the reconstructions with the two operators T (n) : (a) Reconstruction using triple illumination with T (n) C and selection threshold C (n) th = 0.8 Corr n,max . (b) Reconstruction based on energy with T (n) E and selection threshold E (n) th = 0.8 E (n) max . Views at 96 • . White arrows: reconstruction artifacts.

Figure 7 (

 7 b) shows the mean decrease of the ratio between the maximum energy (or correlation) at the current iteration by the maximum energy (or correlation) at the first iteration. The algorithm stops when this ratio is smaller than 0.06.Figures7(a) and 7(b) show as well the cases where the chosen threshold value is too low (E (n) th = 0.5 E (n) max and C

Figure 7 :

 7 Figure 7: (a) Mean decrease in energy of the hologram with time. (b) Temporal variation of the ratio between the mean maximum energy or correlation and its value at the first iteration. The averages of energy, correlation, and time are performed on 256 images. Plots are given for the cases where T (n) considers total energy (red circles) or correlation between the three beams (blue diamonds) as well as for two selection thresholds: 0.8 (empty) and 0.5 (full). Each point corresponds to one iteration of the algorithm.

Figure 8 :

 8 Figure 8: (a) Reconstruction by selection of the sources based on energy with T (n) E and a threshold E (n) th = 0.5 E (n) max . (b) Reconstruction by selection of the sources based on the correlation between the three beams with T (n) C and a threshold C (n) th = 0.5 C (n) max . White arrows: reconstruction artifacts.

Acknowledgments

The authors acknowledge Corinne Fournier and Laurent Daudet for fruitful discussions. Portions of this work were presented at the OSA Digital Holography and 3-D Imaging conference in 2019, paper Th2B.7 [44].

Funding

LabEx Numev (ANR-10-LABX-20: Micro Holo 4D); Centre National de la Recherche Scientifique (CNRS) (Défi: instrumentation aux limites 2018).


Saifeddine Berrayana
email: saifeddine.berrayana@lip6.fr

Etienne Faure
email: etienne.faure@lip6.fr

Daniela Genius
email: daniela.genius@lip6.fr

Frédéric Pétrot
email: frederic.petrot@lip6.fr

Modular On-Chip Multiprocessor for Routing Applications

Simulation platforms for network processing still have difficulties in finding a good compromise between speed and accuracy. This makes it difficult to identify the causes of performance bottlenecks: Are they caused by application, hardware architecture, or by a specificity of the operating system? We propose a simulation methodology for a multiprocessor network processing platform which contains sufficient detail to permit very precise simulation and performance evaluation while staying within reasonable limits of both specification and simulation time. As a case study, we show how a model can be developed for a IPv4 packet routing application, exhibiting the performance and scalability bottlenecks and can thus be used to reason about architectural alternatives.

Introduction

This paper proposes an efficient methodology for the performance tuning of networking applications on network processors by means of describing a flexible platform composed of hardware, software and operating system. Network processors are programmable routers used in specialized telecommunication equipment [START_REF] Shah | Understanding network processors[END_REF]. Among systems on chip they are distinguished by highly parallel hardware: one network processor can contain dozens or even hundreds of microprocessor cores, each executing a massively task parallel application.

In order to arrive at a reasonable performance, the application must be finetuned, and hardware and operating system must be adapted: Our approach is thus typically System-on-Chip. The remainder is structured as follows: First, we present our design methodology. We then show IPv4 as a case study and argue that our platform (hardware and operating system) can be seen as a sharedmemory multiprocessor. Our validation methodology, as well as some performance results are shown before we conclude by outlining current limitations and perspectives.

Description methodology

The network processor modeling platform we took as a reference is ST Microelectronics' STepNP [START_REF] Paulin | STepNP Platform[END_REF], which allows to study network processor architectures on transaction level [START_REF] Cai | Transaction level modelling in system level design[END_REF]. Transaction level modeling reduces simulation time considerably, at the cost of a loss of precision. Our principal aim is to propose a methodology that enables to expose the multiple causes for multiprocessorrelated bottlenecks, hence our need for a fine-grained simulation at register and signal level.

SystemC [START_REF]Open SystemC Initiative: SystemC[END_REF] is a hardware modeling language and simulation kernel that consists of a C++ class library. It permits the instantiation and assignment of models of hardware components and their signals in a hierarchic, building-block manner. Different levels of detail are possible, from coarse grain transaction level to fine grain register transfer level (RTL) [START_REF] Groetker | System Design in SystemC[END_REF].

We opted for a cycle-accurate simulation for which the values present at the connectors of a hardware component are known at every clock cycle. The SO-CLIB component library [START_REF]Projet SOCLIB: Plate-forme de modélisation et de simulation de systèmes integrés sur puce (the SOCLIB project: An integrated systemon-chip modelling and simulation platform[END_REF] provides cycle-accurate, bit-accurate (CABA) simulation models, written in SystemC. Though equivalent to synthesizable models concerning external architecture, these models cannot be directly synthesized -for reasons of protection of intellectual property-and moreover, require much less processing power for their simulation.

The building blocks are easily connected by a unique VCI (Virtual Component Interconnect) interface [START_REF]VSI Alliance: Virtual Component Interface Standard (OCB 2 2.0)[END_REF], which permits to encapsulate arbitrary interconnection networks and to refer to them by a standardized protocol which fixes the number and type of signals required for a communication between components.

Porting an IPv4 Routing Application

A functioning system on chip consists of three parts: The software (the application itself), the hardware and the operating system. Together they build the platform; the challenge is to have them cooperate efficiently. In the next section, we will explain the choices we made to achieve this aim.

Software

IPv4 routing [START_REF] Baker | Requirements for ip version 4 router[END_REF], well-known in the networking community, serves as our reference application. Essentially it consists in taking two or more input streams analyze their headers, and redirect the packets to the appropriate outputs. This application serves as an example in publications on Click [START_REF] Kohler | Click system free software[END_REF][START_REF] Kohler | The Click modular router[END_REF]. Click is a modular router configuration language consisting of two description levels: A simple high-level language to describe the structure of networking applications by composing so-called elements and a library of C++ components containing elementary functionality (such as IP packet header identification, buffering, discarding).

In section 3.3 we will restrict to discussing the application already modified to run on our platform: As can be seen in Figure 1, predefined elements are for example FromDevice and ToDevice, representing entering and exiting packet streams, respectively. Other important elements are Paint and Strip to mark packets for loops and delete a packet header. Within compound elements, several other elements' functionalities are summed up hierarchically. The central functionality of IPv4 lies in such a compound element, LookupIProute, which determines the route taken by an individual packet according to the information in its header. The packets pursue their course along two possible routes: In the direction of ToDevice if the packet is damaged, in the direction of ToSwitch otherwise. A third exit, ToMe, sends ARP requests/responses to the router itself.

Let us now explain the modifications we made for SoC implementation. The first modification concerns I/O. In order to have the application running on one element of a (multiprocessor) network processing platform (one Network Processing Unit, NPU), we have to take into account the fact that its interfaces are asymmetrical (see figure 2). One is the Ethernet link and as such is already treated correctly by Click; the other however is the junction with an on-chip switching matrix. This matrix also uses headers for routing, which have to be constructed and grafted at packet emission and deleted when the packet arrives, respectively.

A consequence of this modification is the addition of these two modules SwitchEncap and FromSwitch which do not appear in the initial list of Click modules. They serve to encapsulate an Ethernet sequence in the header used for internal routing in our chip, and to retrieve an Ethernet sequence after this internal routing, respectively. Moreover we wish to limit complexity. The second modification is due to the structure of the router we use. Our network processor only has two network interfaces connected to the interconnection matrix; likewise IPv4 only has four communication points: Two each for ingress and egress packets.

Thus, it is known that the packets arriving from this interface have already been treated by another block, and are moreover destined for the Ethernet link to the outer world. We can conclude that the verification and routing have already been accomplished by another unit. The only work still to be done is to direct the block to the egress interface. For this reason, Figure 1 shows an additional straight line between fromSwitch and toEth. The only work remaining is to strip the header (delete the first 8 bytes) for the interconnection matrix routing which is expressed in element Strip [START_REF] Baker | Requirements for ip version 4 router[END_REF]. The modified application is shown in Figure 1.

Hardware Architecture Model

The application has been made as independent as possible from the context in which it is used, retaining the essentials of its structure. We look for a modular system architecture where we can easily add, exchange, and regroup hardware models. Such a system is shown in Figure 2. The first step in a modular approach is to define an elementary (SystemC) building block which can be instantiated to the required quantity. This building block is itself a router, and elementary in the sense that it has only two interfaces. This choice enables us to treat platforms with only two interfaces as well as to dimension our platform according to the number of network interfaces required. We have thus a two-level hierarchy of one router composed of four NPUs, interconnected by a switch matrix, each with a certain number of processors inside. From the hardware architecture point of view, two separate parts are to be realized on-chip (we will concentrate on the latter):

an interconnection matrix that links network processors; possible simple but realistic topologies are a full crossbar for a small number of processors, a fat tree for a larger number [START_REF] Andriahantenaina | SPIN: a scalable, packet switched, on-chip micro-network[END_REF] the network processor itself, i e. the unit which does the packet treatment

The hardware architecture we implement was originally inspired by the architecture developed by ST Microelectronics in the context of the STepNP platform. It regroups multiple processors, two network interfaces around a on-chip interconnection network. The platform contains on the one hand calculation units, represented by MIPS R3000 processors with their caches, as well as SRAM memory containing the application code and data; all of these already exist as SOCLIB CABA models. On the other hand, we designed dedicated I/O co-processors to insert Ethernet packets into the system and to extract them, once they are treated. There are two such co-processors for each network interface, one each for ingress and egress, named respectively input-engine and ouput-engine. The architecture is shown in Figure 3. All these elements communicate via an on-chip interconnect carrying a VCI interface. Note that the interconnection network used in the simulations is a virtual model called VGMN (for VCI Generic Micro Network) whose parameters are the number of ports and latency. We are well aware that some precision is lost by using this abstract model instead of a real interconnect model. A first noteworthy advantage of modularity is the fact that the number of co-processors is fixed for each NPU, whatever the number of interfaces that have to be connected to the system. In particular, this avoids to take into account problems due to saturation of interconnect bandwidth in the presence of multiple packet injectors. It will consequently be much easier to extrapolate our performance results to a larger model composed of several instances of the basic model and functioning independently of each other.

The Embedded Software

Once the functionality of our application being fixed, the next step is to port it to the parallel target architecture. Basically, there are two options available: On the one hand, exploiting coarse grain parallelism, which means looking for independent treatments and transform them into separate communicating tasks. On the other hand, duplicating the application in order to obtain identical clones, each treating an IP packet throughout its passage.

The simple kind of routing application we used as benchmark has relatively weak intrinsic parallelism, as one packet is treated by a sequence of functionality corresponding to Click elements. In [START_REF] Chen | Flexible control of parallelism in a multiprocessor PC router[END_REF], the authors propose to decompose one instance IPv4 into threads, with the Fifos as cutting points. As existing routers treat hundreds or thousands of structurally identical tasks in parallel from end to end, we consider it more promising to have all tasks execute the same code on different packets. One RAM contains the global variables, while all other RAMs are allocated to one processor each. However, the accesses to the shared memory are critical, as will show the experimental results.

Operating System Mutek

Once application and hardware fixed, it remains to determine the operating system. Our choice was the Mutek [START_REF] Pétrot | Lightweight implementation of the POSIX threads API for an on-chip mips multiprocessor with VCI interconnect[END_REF] micro kernel. This micro kernel is able to handle multiple tasks running on multiple processors. It provides a C standard library as well as support for POSIX threads.

An important point is that we wish to assign tasks statically to the processors, in order to avoid the cost of task migration between processors, and to avoid migration-related coherence problems. Mutek gives us that possibility.

Once the decisions on kernel and scheduling are made, it remains to determine the number of tasks assigned to on single processing unit. We made performance measures concluding that the time required for context switching outweights time gained by data latency. In the experiments presented here we show one task statically assigned to one processing unit.

Validation and Performance

The hardware and software part of our network processor have to be tested together. The entire platform simulated under SystemC will serve as proof of concept. The validation is made using several Ethernet packets benchmarks injected using the input engine and analyzed by the output engine.

The most relevant performance measure is the throughput in bits per cycle our processor can achieve. Our system is globally synchronous, i e. all of its components share the same clock domain. Thus, we can choose a clock cycle as basic time unit. On the other hand, we have left open the choice of the frequency; clearly, the higher the frequency at which our system works, the better its performance. Our first measure will serve as baseline throughout our experiments: The throughput the system guarantees when one single processor is instantiated, which is the maximal throughput of a single NPU.

When adding processors to our NPU, as can be expected, the performance is not linearly improved. The more adequate question is, how many processors can be added in order to still obtain a profit? Running the simulation during 2 000 000 cycles, and counting the number of 32 bit words arriving at each egress interface allows to neglect the impact of boot time (around 16 000 cycles); also packet(s) that have been treated but not yet been emitted are ignored. All IP packets have 56 bytes, which constitutes the worst case in our application context. This simulation yields in total 689 words of 32 bits each that have been read at egress, a throughput of 0.011 bit/cycle. This rather bad result is due to the non-optimized application. The weak throughput actually corresponds to an average time for treatment of 40 000 cycles for an IP packet. However, keep in mind that our goal is to exhibit performance bottlenecks on a detailed level, for a variety of architectures. Details on the execution times of the individual functions are summed up in table 4 for the case of a packet entering via the Ethernet link and exiting via the interconnection matrix. This corresponds to the longest possible path which a packet can take between two interfaces (around 110 000 cycles). This is nevertheless our reference for multiprocessor experiments. The first consists in simply varying the number of on-chip processing units. Figure 4 illustrates our results, presenting on the x-axis the number of MIPS R3000 processors, on the y-axis the throughput achieved, in bit/cycle. The measurements show that the best throughput is reached for 8 processors, but it is only slightly higher than that for six processors. All units added after the sixth will not yield a significant improvement, they rather lead to performance degradation. The reason for this is that all processors share the same resources, particularly the memory banks where the IP packets are copied to before or after being treated. This contention is only due to the fact that there is only one memory bank for all tasks, not because the data are shared. Thus, we distribute the memory around the interconnection network by replacing the four big memory segments used by the ingress/egress co-processors by as many segments as processors. Such a segment has to fulfill three functions: 1) contain the local data of the thread, 2) allow the ingress processors to write the entering packets there, 3) supply a space for the processor to copy the outgoing packets.

More precisely, each memory bank will be divided into three non overlapping regions. The second region is subdivided into two sections, one for each input engine. As none of the two co-processors knows the behavior of the other, we avoid to make them share variables or address space to not further complicate matters.

This system architecture modification permits a significant performance gain, also for larger numbers of processors. The graph for this second experiment is shown in Figure 5. Here, the performance peak is reached at 27 processors, with a throughput of 0.3 bit/cycle. Here again, a raise in performance to a certain point is followed by an abrupt degradation, caused by contention for the access to shared resources, more precisely, the access to the two segments of shared memory which have remained unchanged and are used by all threads: The code segment and the segment containing the global variables of the system. The presence of contention means that several initiators wish to establish a communication with the target at one time: In consequence, at least one initiator stays blocked waiting for access to the the resource.

The graph in Figure 6 shows the number of cycles due to bottlenecks stemming from accesses to global data in the RAM during an interval of time of 2 000 000 cycles. A closer look reveals that for 28 and more processors, the number of conflict cycles exceeds 1.8 million, which means that the system is paralyzed. Taking a closer look, most of the conflicts happen when processors try to access the RAM containing the operating system data. At the same time, they stay in the scheduler function, which wakes up another thread ready to run when the current thread falls asleep. In our case, there is only one thread per processor, thus we decided to prevent this function from being executed. To do this, the thread should never stop, even if it cannot access the ressource. We changed the lock associated with the shared ressources (input and output engines) from a mutex to a spinlock. From the thread point of view, both have the same behavior: A function call will return only when the resource is available. From the OS point of view, a mutex will suspend the calling thread if the resource is unavailable, whereas a spin lock will continue to demand the resource. More information about mutex and spin locks can be found in [START_REF] Tanenbaum | In: Distributed Operating Systems[END_REF].

New performance results after this software modification are shown in figure 8: The throughput is still growing for more than 27 CPUs, the maximal throughput of 0.56 bit/cycle is reached with 45 processors. Figure 9 compares the three different sets of results we obtained with the different implementations. Obviously, the last change, replacing mutex with spinlocks, does not improve greatly unless we use more than 28 CPUs. This change does not or only marginally improve the throughput per processor, but it allows far more processors to share the same resources without any loss of performances.

Conclusion and Perspectives

By specifying and implementing a simulation model for a network processing platform, we have fulfilled two goals. Firstly, we have determined a detailed and efficient simulation system. The cycle accurate abstraction level is a good compromise between precision and simulation speed. The simulation speed allows us to consider real applications, while the level of detail allows to observe in detail the behavior of our system and helps remedy its insufficiencies. The second goal concerns architecture itself: We have successfully described a network processor core as well as run a multi task application. If performances remain modest, this is mainly due to the software part which will have to undergo profound optimization. Results are on the other hand only slightly impaired by the multiple instantiation of our processing unit. Future work will take several directions: Our next step will be to sum up classes of networking applications by a few "generic" application templates that reflect the typical thread structure -examples for such templates are Quality of Service or Classification applications. Secondly, the use of a more realistic network-on-chip model would improve the precision of our results concerning contention. Thirdly, parts of the multi-threaded micro kernel could be imple-mented in hardware. Finally, our performance studies and the actual implementation have shown the urgent need for debug tools to determine critical resources and time each processor takes for each function. This raises the need for debugging tools for SystemC.

Our platform clearly constitutes a starting point for further experimentation; it is very open in the sense that it allows for a large range of applications from IPv4 routing to traffic analysis and encryption protocols.

Fig. 1 .

 1 Fig. 1. Router Software Configuration

Fig. 2 .Fig. 3 .

 23 Fig. 2. Router with four network interfaces

Fig. 4 .

 4 Fig. 4. Number of processors and resulting throughput in bits/cycle

Fig. 5 .

 5 Fig. 5. Throughput in bits/cycle for the distributed memory implementation as a function of the number of CPUs

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 7. Execution time in cycles required for each element (single-thread, MTU = 56)


Erich Kaltofen
email: kaltofen@math.ncsu.edu

Pascal Koiran
email: pascal.koiran@ens-lyon.fr

Expressing a Fraction of Two Determinants as a Determinant*

Suppose the polynomials f and g in K[x 1 , . . . , x r ] over the field K are determinants of non-singular m×m and n×n matrices, respectively, whose entries are in K∪{x 1 , . . . , x r }. Furthermore, suppose h = f /g is a polynomial in K[x 1 , . . . , x r ]. We construct an s × s matrix C whose entries are in K∪{x 1 , . . . , x r }, such that h = det(C) and s = γ(m+n) 6 , where γ = O(1) if K is an infinite field or if for the finite field K = F q with q elements we have m = O(q), and where γ = (log q m) 1+o(1) if q = o(m). Our construction utilizes the notion of skew circuits by Toda and weakly-skew circuits by Malod and Portier. Our problem was motivated by resultant formulas derived from Chow forms.

Additionally, we show that divisions can be removed from formulas that compute polynomials in the input variables over a sufficiently large field within polynomial formula size growth.

Introduction 1.Motivation

Our investigated problem was motivated by the question of resultant formulas without division. Originally, the resultant of a set of t homogeneous polynomial equations f 1 = • • • = f t = 0 in t variables has been expressed as a GCD of determinants, whose matrices have the coefficients of the polynomials as entries. [START_REF] Macaulay | The Algebraic Theory of Modular Systems[END_REF] gave a formula of a quotient of two such determinants. In special cases, one can remove the division and construct a single determinant that is the resultant [START_REF] Khetan | Determinantal formula for the Chow form of a toric surface[END_REF][START_REF] Eisenbud | Resultants and Chow forms via exterior syzygies[END_REF][START_REF] Khetan | Sylvester A-resultants for bivariate polynomials with planar newton polygons[END_REF]. Those constructions use properties of exact sequences of exterior algebras. However, an algebraic complexity theoretic approach can remove the division in the general case in an entirely different manner. Take Macaulay's formula resultant(f 1 , . . . , f t ) = det(A)/ det(B), where the larger matrix A has dimensions m × m. One converts the determinants to straight-line programs, removes the division by Strassen's [1973] method, or computes their GCD [START_REF] Kaltofen | Greatest common divisors of polynomials given by straight-line programs[END_REF]], parallelizes the straight-line program to O((log m) 2 ) depth [START_REF] Valiant | Fast parallel computation of polynomials using few processors[END_REF][START_REF] Miller | Efficient parallel evaluation of straight-line code and arithmetic circuits[END_REF], converts the resulting division-free straight-line program to a formula of quasi-polynomial size m O(log m) , and finally writes the resultant formula as the projection of a determinant of a matrix C of dimensions k × k where k = m O(log m) [START_REF] Valiant | Completeness classes in algebra[END_REF]]. Note that C with det(C) = det(A)/ det(B) has as entries the coefficients of the original polynomials or constants. Here we shall show that there is a matrix C of dimensions O(m 6 ) whose entries are the coefficients of the original polynomials or constants with det( C) = det(A)/ det(B).

Results and Used Approach

Our main result is the following theorem.

Theorem 1. Let f, g, h ∈ K[x 1 , . . . , x r ] \ {0}, where K is a field, such that f /g = h and f is a projection à la [START_REF] Valiant | Completeness classes in algebra[END_REF] of an m × m determinant and g is a projection of an n × n determinant ( n ≤ m or n > m), meaning that there are matrices A ∈ K[x 1 , . . . , x r ] m×m and B ∈ K[x 1 , . . . , x r ] n×n , whose entries are in K ∪ {x 1 , . . . , x r } with f = det(A) and g = det(B). Then there exists an s×s matrix C whose entries are in K∪{x 1 , . . . , x r } such that h = f /g = det(C) and s is polynomial in m + n, that is, the exact quotient of f and g is a projection of a determinant of polynomial dimension. More precisely,

(i ) if K is infinite or if m = O(|K|) we can take s = O((m + n) 6 ); (ii ) if K is a small finite field, we can take s = O((m + n) 6 • M(log |K| m)), where M(l) = l • (log l) • (loglog l).
We prove our result via the notion of weakly-skew division-free arithmetic circuits by [START_REF] Malod | Characterizing Valiant's algebraic complexity classes[END_REF]. We consider division-free arithmetic circuits (straight-line programs), which are directed acyclic graphs (DAGs) whose nodes have fan-in at most two and which perform addition, subtraction and multiplication. The operands are the values in previous nodes, constant scalars or input variables. The values of designated output nodes are multivariate polynomials in the input variables. The size of the graph is the number of arithmetic operations (sequential complexity) performed. This definition of size is consistent with [START_REF] Valiant | Completeness classes in algebra[END_REF]], for instance. [START_REF] Malod | Characterizing Valiant's algebraic complexity classes[END_REF] work with a slightly different definition of size, in which input nodes (variables or constants) are counted along with arithmetic nodes. To avoid any confusion, we will call this second notion fat size. The fat size of a circuit is therefore equal to the sum of its size and of the number of input nodes. It is bounded by 3 times the size since the number of input nodes is equal to at most twice the number of arithmetic nodes. [START_REF] Toda | Classes of arithmetic circuits capturing the complexity of computing the determinant[END_REF] introduces skew division-free arithmetic circuits, which have the property that at least one of the two operands in each multiplication node is either a scalar constant or an input variable. Toda proves that the determinant polynomial of an m × m matrix can be computed by a skew circuit of size O(m 20 ). In weakly-skew circuits at least one of the two operands to a multiplication node must be computed by a separate circuit. Figure 1 shows an example of a weakly-skew circuits. The separate circuits for the operands of multiplication nodes are marked by dashed boxes.

3 + + + * * x x x x x 1 1 2 1
Figure 1: A weakly-skew circuit [Malod and Portier 2007, Fig. 4]

Valiant [1979] proves that every formula of size s is the projection of an (s + 2) × (s + 2) determinant. In formulas, both operands to all nodes are computed by separate formulas. Valiant's proof can be generalized to show that every polynomial that is computed by a weakly-skew circuit of fat size s is the projection of an (s + 1) × (s + 1) determinant [Malod and Portier 2007, Lemma 6], also using negated input variables. Furthermore, the division-free parallel circuits by [START_REF] Berkowitz | On computing the determinant in small parallel time using a small number of processors[END_REF] for the characteristic polynomial of an m × m matrix can be converted to weakly-skew circuits of size O(m 5 ) [Malod and Portier 2007, Proposition 5]. In Section 3 we apply Strassen's [1973] technique for elimination of divisions directly to the characteristic polynomials derived from A and B in Theorem 1 above, thus obtaining a division-free weakly-skew circuit for h = det(A)/ det(B), which then is the projection of a determinant. Note that the technique in [START_REF] Canny | Generalized characteristic polynomials[END_REF]], which is not directly applicable, assumes that the diagonals of A and B hold a single separate variable, implying also n < m. Along the way we show in Section 2 that weakly-skew circuits can be simulated by skew circuits with an increase in size by a constant factor only; and that Toda's O(m 20 ) bound can be reduced to O(m 5 ).

Since weakly-skew circuits are projections of determinants, we have the following corollary to Theorem 1 and our transformation results to skew circuits.

Corollary 1. Let f, g, h ∈ K[x 1 , . . . , x r ] \ {0}
, where K is a field with deg(g) = O(|K|) such that f /g = h and f and g are computed by a weakly-skew circuit of size s with inputs x 1 , . . . , x r . Then h can be computed by a skew circuit of size O(s 6 ). Suppose in the above corollary that f and g are computed by division-free formulas of size ≤ s. Then h = f /g is a projection of a determinant of a matrix of dimension O(s 6 ). In Section 4 we show that there exists a division-free formula of size (s 5 ) that computes h = f /g. As a consequence of the latter result, we can show that a polynomial in K[x 1 , . . . , x r ] that is computed by a formula of size s with additions, subtractions, multiplications and divisions, where K is a sufficiently large field and x 1 , . . . , x r are input variables, can be computed by a division-free formula of size s O(1) (see Theorem 3 on page 10).

From Weakly-skew to Skew Circuits

In this section we show that weakly-skew circuits can be efficiently simulated by skew circuits. This fact will be used in the proof of Theorem 1, part (ii) given in Section 3.2.

Two algorithms for converting weakly-skew circuits into skew circuits are already described by [START_REF] Malod | Characterizing Valiant's algebraic complexity classes[END_REF]. The polynomial computed by a weakly-skew circuit of fat size m can be converted into a determinant of size m + 1 by [Malod and Portier 2007, Lemma 6]. One can then apply Toda's algorithm, which evaluates a determinant of an m×m matrix by a skew circuit of size O(m 20 ). It is observed in [Malod and Portier 2007, Section 5.2] that the role of the determinant in this first algorithm can be played by the polynomial F m = Trace(X m ). The method described below is even simpler, and more efficient. As a byproduct, we obtain an improvement to O(m 5 ) of Toda's original O(m 20 ) bound.

Remark 1. In the determinant constructed in [Malod and Portier 2007, Lemma 6], all input variables are negated. This is not a problem for converting a weakly-skew circuit into a skew circuit with the algorithm described in the paragraph above. However, the occurence of negated input variables makes their construction unsuitable for the proof of Theorem 1. To circumvent this difficulty, one can modify slightly Malod and Portier's condition to obtain a matrix of size m + 1 without negated variables. This can be done in two steps:

1. In their proof of Lemma 6 Malod and Portier first construct a matrix B of size m such that det B = -f , where f is the polynomial computed by a weakly-skew circuit of size m. In this matrix all variables are negated (and all diagonal entries are equal to 1 except the first, which is 0). Now let C = -B: we have det C = ±f , and there are no negated variables in C.

2. If m is odd we have det C = -f . As Malod and Portier, we can add one last row and one last column full of 0's (except for an entry equal to -1 in the bottom right corner) to obtain a matrix of size m + 1 whose determinant is equal to f .

We shall work with acyclic edge-weighted directed graphs. We recall that the weight of a path in such a graph is defined as the product of the weights of the edges appearing in the path. If s and t are two vertices of G, the weight of (s, t) in G is defined as the sum of the weights of all paths from s to t.

Lemma 1. Let W be a weakly-skew circuit of fat size m. There exists an acyclic directed graph G, with two distinguished vertices s and t such that:

(i ) The weight of (s, t) in G is the polynomial computed by W , and G is of size at most m + 1.

(ii ) Every vertex in G other than s has either a single incoming edge, of weight equal to an input of W or to the constant 2, or two incoming edges, each of weight 1. This is essentially Lemma 5 of [START_REF] Malod | Characterizing Valiant's algebraic complexity classes[END_REF]. In that lemma the authors prove the existence of a graph G satisfying (i). An inspection of their proof shows that the graph that they construct also satisfies the second property.

Proposition 1. Let W be a weakly-skew circuit of fat size m. There exists a skew circuit W ′ which is equivalent to W (i.e., computes the same polynomial), has the same number of input nodes as W , and has at most m arithmetic nodes (W ′ is therefore of fat size at most 2m).

The algorithmic idea behind this result is quite simple: for each vertex v in the graph G of Lemma 1 we compute the weight ω(v) of the pair (s, v). If v has a single incoming edge of weight x connecting s to v then of course ω(v) = x. If v has a single incoming edge of weight x connecting a vertex v ′ = s to v, we can apply the formula ω

(v) = x × ω(v ′ ) if x is an input of W . If x is the constant 2, we apply the formula ω(v) = ω(v ′ ) + ω(v ′ ) (one could of course use the the formula ω(v) = 2 × ω(v ′
) instead, at the cost of introducing one additional constant input in W ′ ). Finally, if v has two incoming edges of weight 1, connecting the vertices v 1 and v 2 to v, we have ω(v) = ω(v 1 ) + ω(v 2 ). The resulting circuit satisfies the requirements of Proposition 1.

Corollary 2. The determinant of a m × m matrix can be computed by a skew circuit of size O(m 5 ). This is a significant improvement over the O(m 20 ) bound given by [START_REF] Toda | Classes of arithmetic circuits capturing the complexity of computing the determinant[END_REF]. Corollary 2 is an immediate consequence of Proposition 1 since, as pointed out in section 1, the determinant of a m × m matrix can be computed by a weakly-skew circuit of size O(m 5 ).

Skew circuits will be useful in Section 3.2 due to the following proposition and the subsequent remarks.

Proposition 2. Let K be a field and let L = K(θ) be an algebraic extension of K of degree d = [L : K]. Let f ∈ L[x 1 , . . . , x r ] be a polynomial computed by a skew circuit W L of size m, with input nodes labeled by variables from {x 1 , . . . , x r } or constants from L. Let us expand f according to the powers of θ: one can write f = d-1 j=0 θ j f j , where f j ∈ K[x 1 , . . . , x r ]. There exists a skew circuit W K of size O(d 2 m) with at most d output nodes which computes simultaneously all the polynomials f j . Moreover, W K uses only constants from K.

Proof. Let α be a node of W L computing a polynomial f α ∈ L[x 1 , . . . , x r ]. We use a standard technique: in W K , f α will be represented by d nodes computing polynomials f 0,α , . . . ,

f d-1,α ∈ K[x 1 , . . . , x r ] such that f α = d-1 j=0 θ j f j,α .
If α is an input node labeled by some variable x i we can take (f 0,α , . . . , f d-1,α ) = (x i , 0, . . . , 0); if α is labeled by the constant d-1 j=0 a j θ j we can take (f 0,α , . . . , f d-1,α ) = (a 0 , a 1 , . . . , a d-1 ). Assume now that α is an addition node with inputs coming from nodes β and γ. In this case we simply perform componentwise additions since f j,α = f j,β + f j,γ .

Finally, the case where α is a multiplication node can be split in two subcases since W L is skew: multiplication by a variable x i , or multiplication by a constant from L. In the first subcase, assume that α multiplies the output of node β by x i . We have f j,α = x i f j,β . Observe that the d resulting multiplications are skew. In the second subcase, assume that α multiplies the output of node β by a constant η ∈ L. Multiplication by a constant is a Klinear operation. The tuple (f 0,α , . . . , f d-1,α ) can therefore be obtained from (f 0,β , . . . , f d-1,β ) by multiplication by an appropriate d × d matrix A η with entries in K. The corresponding matrix-vector product can be computed in O(d 2 ) arithmetic operations, and once again the resulting multiplications are all skew.

Remark 2. We will apply this result in Section 3.2 in a situation where we know that the output of W L lies in fact in K[x 1 , . . . , x r ]. In this case, the subcircuit associated to the first output node of W K computes the same polynomial as W L .

Remark 3. Proposition 2 also applies to weakly-skew rather than skew circuits. First one converts the weakly-skew circuit into a skew circuit using Proposition 1. One can then apply Proposition 2 to the skew circuit.

Elimination of Divisions

Large Coefficient Fields

Our symbolic determinants live in a multivariate polynomial domain K[x 1 , . . . , x r ], where K is a sufficiently large field. Consider we have non-singular matrices A(x 1 , . . . , x r ) ∈ K[x 1 , . . . , x r ] m×m and B(x 1 , . . . , x r ) ∈ K[x 1 , . . . , x r ] n×n , whose entries are either variables or constants, i.e., ∀i, j, k, l with 1

≤ i, j ≤ m, 1 ≤ k, l ≤ n : (A) i,j , (B) k,l ∈ K ∪ {x 1 , . . . , x r }.
Here (M ) i,j denotes the element in row i and column j in the matrix M .

We suppose now that det(A)/ det(B) ∈ K[x 1 , . . . , x r ], that is the polynomial division by det(B) is exact. We construct a division-free weakly-skew arithmetic circuit W of size O((m + n) 6 ), i.e, polynomial in the dimensions of A and B, that computes the polynomial det(A)/ det(B).

The construction follows Strassen's [1973], using Berkowitz's [1984] ‡ and [START_REF] Chistov | Fast parallel calculation of the rank of matrices over a field of arbitrary characteristic[END_REF] weakly-skew arithmetic circuits for the characteristic polynomial. Let u 1 , . . . , u r ∈ K be such that both U A = A(u 1 , . . . , u r ) and U B = B(u 1 , . . . , u r ) are non-singular. Such values always exist if |K| > m ≥ deg(det(A)) [START_REF] Schwartz | Fast probabilistic algorithms for verification of polynomial identities[END_REF][START_REF] Zippel | Probabilistic algorithms for sparse polynomials[END_REF]]. ‡ Sasaki and Murao [1982] compute the characteristic polynomial of an n × n matrix with entries in a commutative ring in n ω+1+o(1) ring operations, and Berkowitz [1984, Section 4] conjectures that there exists a division-free arithmetic circuits of size O(n ω ) for the characteristic polynomial, given that n × n matrices can be multiplied with O(n ω ) operations. We can set ω = 2.375477 [START_REF] Coppersmith | Matrix multiplication via arithmetic progressions[END_REF]. The Sasaki&Murao/Berkowitz problem remains open. The best division-free complexity for the characteristic polynomial is O(n 2.697263 ) [START_REF] Kaltofen | On the complexity of computing determinants[END_REF].

Now consider det(U

B -λ(U B -B)) = det(U B ) • det(I n -U -1 B λ(U B -B) λ B ) ∈ K[x 1 , . . . , x r ][λ], (1)
where I n denotes an n-dimensional identity matrix. Now the coefficient of λ i in (1) is the homogeneous part, hom T (det(B), i), of total degree i of det(B) represented in the term basis

T = {(u 1 -x 1 ) d 1 • • • (u r -x r ) dr | d j ≥ 0}, (2) namely det
(B) = n i=0 d 1 +•••+dr=i c d 1 ,...,dr (u 1 -x 1 ) d 1 • • • (u r -x r ) dr hom T (det(B), i)
where c d 1 ,...,dr ∈ K.

Note that evaluation (1) at λ = 1 gives det(B). We compute 1/ det(U B -λ(U B -B)) as a truncated power series in

K(u 1 -x 1 , . . . , u r -x r )[[λ]].
Because the constant coefficient of (1) as a polynomial in λ is in K, the coefficients of λ i in the power series for the reciprocal are homogeneous polynomials of degree i in the basis (2). We present a weakly-skew circuit for the coefficients of λ i . Let M 1...l,1...l denote the top left l × l principal submatrix of a matrix M . Chistov's algorithm is based on the identities

1 det(I n -λ B) = n j=1 det(I j-1 -λM 1...j-1,1...j-1 ) det(I j -λM 1...j,1...j ) = n l=1 (I l -λ B1...l,1...l ) -1 l, l = n l=1 ∞ k=0 λ k Bk 1...l,1...l l, l ≡ n l=1 m k=0 ( Bk 1...l,1...l ) l, l λ k (mod λ m+1 ) (3) ≡ 1 + q 1 λ + • • • + q m λ m (mod λ m+1 ). (4)
We use weakly-skew circuits to compute the coefficients q k . Each ( Bk

1...l,1...l ) l, l in (3) is com- puted as B1...l,1...l • ( B1...l,1...l • (• • • ( B1...l,1...l    0 . . . 0 1    ) . . .))
l by weakly-skew circuits of size O(mn 2 ). The weakly-skew circuits for carrying out the modular product (4) require no more than m + 1 copies of the circuits for each coefficient ( Bk 1...l,1...l ) l, l . Thus a weakly-skew circuit W 1 of size O(m 2 n 3 ) computes all q k in (4). Malod and Portier [2007, Proposition 5] compute the coefficients of

1 + p 1 λ + • • • + p m λ m = det(I m -U -1 A λ(U A -A))
via a weakly-skew circuit W 2 of size O(m 5 ) by Berkowitz's algorithm. Alternatively and less efficiently, one could as above compute

1 det(I m -U -1 A λ(U A -A)) ≡ 1 + p1 λ + • • • + pm λ m (mod λ m+1 ),
and compute the truncated power series of the reciprocal as

1 1 + p1 λ + • • • + pm λ m ≡ m l=0 (-p 1 λ -• • • -pm λ m ) l (mod λ m+1 ) (5) ≡ 1 + p 1 λ + • • • + p m λ m (mod λ m+1 ),
again repeating the circuits which compute the coefficients pi in the truncated polynomial powers.

A weakly-skew circuit then carries out the multiplication

det(U A -λ(U A -A)) det(U B -λ(U B -B)) ≡ det(U A ) det(U B ) (1 + p 1 λ + • • • + p m λ m ) • (1 + q 1 λ + • • • + q m λ m ) (mod λ m+1 ) (6) ≡ ̺ 0 + ̺ 1 λ + • • • + ̺ m λ m (mod λ m+1 ).
Here and before, the truncation could be performed at deg(det(A)/ det(B)) + 1 ≤ m + 1. Again, we need no more than m + 1 copies of the circuit W 1 for the coefficients q k or of W 2 for p k . Because h = det(A)/ det(B) has total degree ≤ det(A) ≤ m and ρ i = hom T (h, i), we can compute det(A)

det(B) = ̺ 0 + ̺ 1 + • • • + ̺ m
by a division-free weakly-skew arithmetic circuit of size O(min{m 3 n 3 + m 5 , m 2 n 3 + m 6 }).

Remark 5. The above techniques also yield a single determinant of a matrix of polynomiallysized dimensions for a fraction of products of determinants ( i A [i] ) ( j B [j] ) ∈ K[x 1 , . . . , x r ], with and without using block diagonal matrices. Such fractions occur when computing the resultant via Koszul complexes § .

Small Coefficient Fields

When the coefficient field K has few elements, divisions by zero may occur at all values u 1 , . . . , u r ∈ K and there is no no-singular U B in Section 3. We can nonetheless obtain the bound of Theorem 1.(ii) using a field extension. Indeed, from Section 3 we know that the quotient h = f /g can be computed by a polynomial size weakly-skew circuit W L with constants from L if L is an extension of K with § Ágnes Szántó has pointed this application out to us.

at least m + 1 elements. We can therefore work with an extension of degree O(log |K| (m)). Since h is actually a polynomial with coefficients in K, by applying Proposition 2 and the two subsequent remarks to W L we obtain a skew circuit W K which computes h using constants from K only. Finally, W K can be transformed into a determinant by [Malod and Portier 2007, Lemma 6] as in the proof of Theorem 1(i).

Formulas With Divisions

We shall prove the following theorem.

Theorem 2. Let f, g, h ∈ K[x 1 , . . . , x r ] \ {0}, where K is a field with |K| > deg(f ) + deg(g) × deg(h)
, such that f /g = h and f is computed by a division-free formula of size s f and g is computed by a division-free formula of size s g with inputs x 1 , . . . , x r . Then h can be computed by division-free formula of size O((s f + s g ) 5 ).

The proof uses interpolation (cf. [Kaltofen 1988, Section 5]). Denote by the total degrees δ f = deg(f ), δ g = deg(g) and δ h = deg(h). Note that δ f ≤ s f and δ g ≤ s g . As in Section 3, let u 1 , . . . , u r ∈ K such that g 0 = g(u 1 , . . . , u r ) = 0. Again as in ( 5) and ( 6) using the term basis (2), we compute

δ h i=0 hom T (h, i)λ i ≡ δ f i=0 hom T (f, i)λ i δg i=0 hom T (g, i)λ i (mod λ δ h +1 ) ≡ 1 g 0 δ f i=0 hom T (f, i)λ i • δ h l=0 - 1 g 0 l δg i=1 hom T (g, i)λ i l H(x 1 , . . . , x r , λ) (mod λ δ h +1 ). (7)
We compute the polynomial H in (7) of degree in λ of no more than D = δ f + δ g δ h by interpolation at λ = v 0 , . . ., v D ∈ K as an exact polynomial before truncating modulo λ δ h +1 . Note that for all j we have

δ f i=0 hom T (f, i)v i j = f (u 1 -v j (u 1 -x 1 ), . . . , u r -v j (u 1 -x r )) and δg i=1 hom T (g, i)v i j = g(u 1 -v j (u 1 -x 1 ), . . . , u r -v j (u 1 -x r )) -g 0 , (8)
and therefore can obtain the values in the inputs x 1 , . . . , x r by formulas. We have formulas for each H(x 1 , . . . , x r , v j ) by repeating the formulas (8) no more than δ h (δ h + 1)/2 many times.

Interpolation is a matrix times vector product and again is done by repeating the formulas for H(x 1 , . . . , x r , v j ) no more than D + 1 times. Finally, we add the thus obtained first δ h + 1 coefficients in λ of H. Note that all divisions are by scalars independent on x 1 , . . . , x r .

Theorem 2 together with the well-known parallel circuits for formula evaluation allows the removal of divisions in formulas altogether. When there are division nodes in formulas with inputs x 1 , . . . , x r , it is assumed that all rational functions in K(x 1 , . . . , x r ) by which is divided are non-zero. Formulas with a divisions by a generic 0 are naturally excluded. For certain values in K for the inputs x 1 , . . . , x r a zero division can occur.

Theorem 3. There exists a real constant γ > 0 with the following properties. Let h ∈ K[x 1 , . . ., x r ] be computed by a formula (with divisions) of size s with inputs x 1 , . . . , x r . Assume K is a field with |K| > s γ . Then h can be computed by division-free formula of size O(s γ ).

In the following proof γ 1 , γ 2 and γ 3 are fixed positive real constants. The proof observes that h is computed by a circuit V 1 of fan-in at most 2 with divisions of depth ≤ γ 1 log(s) ([ [START_REF] Kosaraju | Optimal parallel evaluation of tree-structured computations by raking[END_REF] and the references there). By computing unreduced numerator and denominator polynomials for each node separately, we have two division-free circuits V 2 and V 3 of depth ≤ γ 2 log(s) that compute polynomials f and g such that f /g = h. We can convert V 2 and V 3 into division-free formulas of depth ≤ γ 2 log(s), hence of size < 2 γ 2 log(s)+1 = O(s γ 3 ), which also bounds the degrees of f and g. Applying Theorem 2 to both formulas yields Theorem 3.

Remark 4 .

 4 The size of W K in Proposition 2 can be reduced to O(m d (log d) loglog d) by fast polynomial multiplication algorithms [Cantor and Kaltofen 1991] and fast division with remainder algorithms [von zur Gathen and Gerhard 1999, Section 9.1].

* This material is based on work supported in part by the National Science Foundation under Grant CCF-0514585 (Kaltofen).