Alban Fichet

Romain Pacanowski

Alexander Wilkie

Marc Stamminger

An OpenEXR Layout for Spectral Images

Bi-spectral image

Examples of using the proposed OpenEXR layout to store spectral and bi-spectral information along with the polarization state. Left: For any given pixel of a radiance image, we can extract the full spectral data (blue box). We also store the corresponding Stokes components per wavelength, shown here averaged for all stored wavelengths in a specific area of the image (red box). Right: A bi-spectral image can store a full re-radiation matrix for each pixel. This is obviously only meaningful as input data for a reflective texture. The sample points are two locations with different bi-spectral reflectance, fluorescent pink (blue box) and fluorescent yellow (red box). The corresponding RGB image uses a D65 illuminant. These two OpenEXR spectral images are provided in the supplemental materials.

Introduction

In recent years, spectral rendering has become increasingly popular, e.g., for high end visual effects work [START_REF] Fascione | Manuka: A batch-shading architecture for spectral path tracing in movie production[END_REF]]. The main reason for its increasing adoption is that it allows for a far more physically plausible handling of various important vi-1 ISSN 2331-7418 sual effects such as skin appearance, at least when compared to tristimulus systems. It also allows proper integration of still rarely found features such as dispersion, polarization, and fluorescence. Spectral renderer output is not tied to any particular color space, whereas tristimulus images are intrinsically only valid for a given combination of white point and RGB primaries. However, the output of spectral rendering systems such as Manuka [START_REF] Fascione | Manuka: A batch-shading architecture for spectral path tracing in movie production[END_REF]] is usually not stored in spectral form: the main reasons why this is rarely done are spectral aliasing issues, storage considerations, and the fact that, for most workflows, wide-gamut RGB data in floating point format is sufficient from the point of rendering onward. However, there are applications in both research and production environments where storing the actual spectral output of a renderer is useful, so we still need a standardized way of doing this. And there is one additional application area where the ability to use spectral images is arguably far more important anyway: as input for the shading subsystem of a spectral renderer. If accurate matches between measured surface appearance and renderings are desired, the ability to handle (bi-)spectral reflectance textures is essential.

Many RGB image formats exist, but-somewhat surprisingly-the same is not true for spectral images. Each pipeline uses either its own proprietary flavor of image format (e.g., ENVI hyperspectral [Harris Geospatial Solutions 2020]) or a collection of monochromatic images stored in a set of RGB image formats (e.g., META-COW [START_REF] Fairchild | METACOW: A public-domain, highresolution, fully-digital, noise-free, metameric, extended-dynamic-range, spectral test target for imaging system analysis and simulation[END_REF]). This makes the exchange, support in software, and long-term storage of spectral images difficult.

We propose a spectral image format based on the OpenEXR framework: the result is suitable not only for rendering but also for various applications such as acquisition, and as a storage format for (bi-)spectral reflectance. In addition, we retain the user friendliness of an RGB image format allowing a quick preview of the spectral image in commonly used graphics software.

Related Work

To store spectral images, a simple technique is to provide a collection of standard RGB image files, where each image corresponds to a specific spectral band [START_REF] Yasuma | Generalized assorted pixel camera: Postcapture control of resolution, dynamic range, and spectrum[END_REF]]. The key advantage of this method is that it allows uncomplicated viewing of all spectral channels with standard software. On the other hand, handling multiple separate files is not convenient, and reconstruction of color images from the spectral data is tedious.

Specialized formats exist for storing spectral data in specialized fields. For example, the ENVI format [Harris Geospatial Solutions 2020] targets satellite imagery, and separate data blocks are stored for each spectral channel. A header file, in plain text format, describes the information stored in the images, and a binary file is used as a dump of the image data. This reduces the amount of separate file types, but still requires two separate files in order to read the image. Moreover, there is no standardized way to store polarization and bi-spectral information, as especially the latter never occurs in satellite imaging contexts.

In the rendering context, the ARTRAW format was introduced along with the Advanced Rendering Toolkit (ART) renderer [Advanced Rendering Toolkit 2018]. It supports polarized emissive images, and the spectral sampling options are those of the renderer itself (8, 11, 18 and 46 spectral samples, respectively). The header and data are stored in a single file, in a mixture of plain text (header) and binary (data) formats. This format is mostly used only by this specific renderer and cannot be displayed without specialized software: it also intentionally eschews all attempts at reducing file size, as it was only ever intended as a lossless reference data storage for raw rendering output. As spectral images tend to be quite large, any reduction in file size would be very welcome; so, adopting ARTRAW (which is not very flexible in terms of its internal structure anyway) as a standard would be decidedly suboptimal.

Finally, there is OpenEXR [START_REF] Kainz | The OpenEXR image file format[END_REF]], the de facto standard for storing HDR color space images. The official OpenEXR specifications do not cover the use and description of spectral data, but allow very flexible storage of multiple layers of binary data. Some software packages already use layers of the OpenEXR format to store spectral data [Eclat Digital 2019], but, to the best of our knowledge, the specifications used there are limited. These formats generally lack support of polarization, bi-spectral, and capture metadata information with a unified interface.

Scope of Our Spectral Image Representation

In this section, we describe the range of information we intend to store with our spectral representation.

Reflective and Emissive Spectra

First, we distinguish between emissive and reflective spectra. Both are used in distinct places within a rendering system, and do not share the same physical units.

An emissive spectrum L e (λ) describes radiant energy within a scene, such as the radiance emitted by light sources. The image computed by a path tracer is a collection of emissive spectral values: it represents the spectral radiance reaching the sensor. The pixels in an environment map are a specialized case of light source and also contain emissive spectral data. L e (λ) is expressed in W • m -2 • sr -1 • nm -1 in radiometric units. 1 On the other hand, a reflective spectrum describes the attenuation of energy af-ter reflection, scattering, or transmission/volume absorption events. For instance, an albedo texture image describes how a surface attenuates incoming radiant energy upon reflection. The albedo is unitless and is normally used in the context of a BRDF model (in sr -1 ) that further describes the directional dependency of surface reflectance:

ρ λ (ω o ) = θ i φ i brdf λ (ω o , θ i , φ i ) cos θ i sin θ i dθ i dφ i [sr -1 ].
(1)

Polarization

Although it is not directly visible to human eyes, 2 the polarization state of electromagnetic radiation can play an important role in light transport for some scenes. Specular inter-reflections are mis-predicted to varying degrees by non-polarizing rendering computations, as well as the reflections of blue sky radiance off specular surfaces.

For various purposes related to rendering, such as being able to store a polarized environment map of sky dome radiation and to inspect the polarization state of reflective highlights in a scene, we want a spectral image format to optionally also store polarization information.

The mathematical formalism of Stokes vectors and Mueller matrices allows the representation of fully, partially, or non-polarized light in a compact and intuitive manner that blends in well with non-polarized spectral rendering. Stokes vectors describe light polarization,

S =      S 0 S 1 S 2 S 3      , (2)
and Mueller matrices describe the effect of a scattering event on light polarization: it encodes the transform of a Stokes vector [START_REF] Goldstein | Polarized Light[END_REF][START_REF] Wilkie | How to write a polarisation ray tracer[END_REF].

In order to achieve the desired functionality for emissive images, we store the following information, per pixel and per wavelength:

• all Stokes components S 0 , S 1 , S 2 , S 3 if the image contains polarization information, or

• only the first Stokes component S 0 if the image does not contain polarization information.

To use a consistent notation, we always refer to emissive spectral radiance in images as Stokes component S 0 , even in non-polarized images.

Bi-spectral Reflectance

Even though this behavior is not yet commonly simulated in current rendering software, a significant number of manmade and natural materials exhibit wavelength shifting behavior when reflecting or scattering light. Apart from obviously fluorescent artificial materials such as Day-Glo warning paint, examples that show at least some degree of such behavior include wood (lignin is fluorescent), green plants (chlorophyll), and paper and textiles (mostly due to washing powders usually containing optical brighteners that make white cloth look "whiter" by adding a blue fluorescent tinge to it). Similar to storing "normal" spectral albedo textures that are then used in the context of a BRDF model, complex appearance modeling requires the ability to store extended albedo textures that also describe wavelength shifting. As it can be seen in Figure 2, one needs to store a two-dimensional spectrum-the re-radiation matrixfor each pixel in this case: for a given incoming wavelength λ i , re-emission can occur at various intensities for a wide range of possible output wavelengths λ o .

Capture Information

For measured spectral images, we also need to store the capture device characteristics as metadata. The Cornell Box website [Cornell University Program of Computer Graphics 1998] serves as an example of the sort of data we need to retain: as reference data, they provide a set of monochromatic images of the real Cornell Box from 400 nm to 700 nm with a sample taken every 50 nm.

In the Cornell dataset, each image was taken with a monochrome camera with a specific narrow-band filter placed in front of the lens: detailed information about the measurement setup is provided along with the raw data on the website. If one wants to be able to interpret the contents of such a measured image in a stand-alone form, directly from the image, this information should be included in the image header. There are two optical elements that alter the energy received by the sensor: First, the camera lens has a characteristic transmission spectrum. Second, each measurement filter has a transmission spectrum. In addition to this, the camera sensor response is usually not linear and also needs to be stored. If all this information is retained in the image header, a captured spectral image can be used for long-term archival storage without losing the ability to properly interpret its contents.

Using OpenEXR for Spectral Data

OpenEXR supports the storage of multiple layers within a single file, but there are some constraints on using this functionality. In this section, we explain our choices made in order to meet our requirements for storing spectral data without violating the OpenEXR format specification.

Layers and Channels Organization

The official OpenEXR specification [START_REF] Kainz | Reading and writing OpenEXR image files with the IlmImf Library[END_REF]] states that an OpenEXR image can consist of different layers and that each layer can contain an arbitrary number of channels. In other words, a layer is a group of channels, and from an implementation point of view, layers and channels are indistinguishable. They are merely separated by a naming convention that is enforced by convenient functions in the official OpenEXR library. A layer name is followed by a dot character. Layers can be nested and the tailing string names a channel.

The following outlines the layer naming conventions that we propose. length/frequency values are always positive values, and an optional integer power of ten exponent also can be provided with E or e followed by an optional sign and an integer value. This value is followed by an optional unit multiplier and by the used unit. These are either Hertz (Hz) or meters (m). Table 1 enumerates the different fields that describe a channel name. The channel naming convention is also described in Appendix A.2.

For instance, the channel "S0.560,5nm" corresponds to an S 0 component at 560.5 nm.

A regex matching this naming convention can be written as follows: )" units = "(Y|Z|E|P|T|G|M|k|h|da|d|c|m|u|n|p|f|a|z|y)?(m|Hz)" regex = "ˆ" + layer + '\.' + value + units + "$"

Bi-spectral Data

To store a re-radiation matrix for each pixel of a reflective image, we can also use the OpenEXR layering system, although the resulting image is of considerable size. The main diagonal, which corresponds to the case where there is no re-radiation, follows the same convention as described in Section 4.1.2. Nested layers store the re-radiation part named with the re-radiation wavelength or frequency using the same comma separated format followed by the units. For instance, re-radiation from 560 nm to 600 nm would be stored in a channel named "T.560nm.600nm."

Nested Hierarchy

We also define a layer convention for child layers, so that it is possible to use a spectral layout within an existing layer hierarchy. For instance, a stereo image with 550 nm layers will have its "default" (left) spectral layer stored in "S0.550nm" and its right view in "right.S0.550nm."

Vol. 10, No. 3, 2021 http://jcgt.org

Metadata

We allow the storage of additional metadata to retain capture information using OpenEXR custom attributes in the image header. OpenEXR documentation recommends using types provided from ImfStandardAtrributes.h. A summary of all metadata we support is provided in Appendix A.1.

Layout Version

To allow future evolution of the proposed spectral layout, we add an attribute that specifies the format version. If the layout is exactly as described in this article, the value shall be set to "1.0" (see Table 2).

Description Attribute Name Value Version spectralLayoutVersion "1.0"

Radiometric Units

When the image contains emissive spectral layers, the header has to specify the used radiometric units (see Table 3). For instance, the result of a rendering or a hyperspectral image shall be in

W • m -2 • sr -1 .

Description

Attribute Name Value Emissive layers radiometric units emissiveUnits "W," "W.mˆ-2," "W.srˆ-1," or "W.mˆ-2.srˆ-1"

Table 3. Emissive unit information in the OpenEXR header: this is only mandatory for emissive images.

Polarization

When using emissive polarized images, we need to specify the orientation convention used (see Table 4). Figure 3 illustrates how the same signal can be described as left circular polarized or right circular polarized depending on the convention used.

Description

Attribute Name Value Polarization frame polarisationHandedness "left" or "right" And the special cases for the first and the last bins are

B 1 (λ -) = B 1 - B 2 -B 1 2 , (5)
B N (λ + ) = B N + B N -B N -1 2 . (6)
Figure 4 illustrates this default filtering behavior with regular and varying sampling rate.

Spectrum Type Metadata

To store spectrum types, we use an Imf::StringAttribute packing spectrum in a string. Each pair of wavelength or frequency and value are separated by a colon and followed by a semicolon. We enforce specification of the units of the wavelength or frequency and allow a unit multiplier. Contrary to channel naming, we do not have a constraint on the character set, so we use dot decimals.

A regex matching this naming convention can be written as follows:

value = "(\d * \.?\d * ([eE][-+]?\d+)?)" units = "(Y|Z|E|P|T|G|M|k|h|da|d|c|m|u|n|p|f|a|z|y)?(m|Hz)" regex = value + units + ":" + value + ";"

For instance, a triangular filter response from 620.5 nm to 679.5 nm with its maximum at 650 nm is described as follows:

Data Compression

The OpenEXR format offers a wide variety of lossless and lossy compression algorithms, which is interesting for spectral and bi-spectral data that have a significant memory footprint. An overview of the various compression formats available are discussed in the technical introduction to OpenEXR [START_REF] Kainz | Technical Introduction to OpenEXR[END_REF]]: choosing any particular one of them depends on the application one needs spectral image storage for (that is, whether lossy compression is acceptable, or not). When in doubt, lossless compression should be chosen.

Enforcing Backward Compatibility

OpenEXR is widely supported in graphic packages (GIMP, Photoshop, etc.) in its RGB variant. Although these software packages only support normal color spaces (RGB, CMYK, etc.), we can provide additional layers in spectral OpenEXR files that allow opening and previewing an RGB version of the spectral image (see Figure 5). This is simply done by adding R, G, and B channels to the layer set. This behavior, which is not necessary to have a valid spectral OpenEXR file, is strongly encouraged in the OpenEXR official documentation page [START_REF] Hillman | Storing multi-view images in OpenEXR files[END_REF].

Multiple spectral-to-RGB conversions exist. In this section we explain how we keep track of the conversion pipeline.

Vol. 10, No. 3,2021 http://jcgt.org

Proposed Conversion Pipeline

Emissive Images

For emissive images, we use the S0 layer and convert the spectrum to RGB using the following pipeline:

1. Convert emissive spectrum L e (λ) to CIE-XYZ 2 • 1931.

2. Convert XYZ to linear sRGB D65.

Steps 1 and 2 correspond to the following equation:

   R G B    = M XYZ→RGB •    Λ L e (λ) • x(λ)dλ Λ L e (λ) • ȳ(λ)dλ Λ L e (λ) • z(λ)dλ    , (7)
where

• M XYZ→RGB =    3.2404542 -1.5371385 -0.4985314 -0.9692660 1.8760108 0.0415560 0.0556434 -0.2040259 1.0572252    ,
• L e (λ) is the emissive spectrum,

• x(λ), ȳ(λ), z(λ) are the sensitivity curves for the CIE-XYZ 2 • 1931 color space.

Reflective Images

For a reflective image, we use the T layer and convert the spectrum to RGB using the following pipeline:

1. Multiply the spectrum S(λ) with CIE D65 illuminant L D65 (λ). Spectral case. The conversion is defined as follows:

   R G B    = 1 Y D65 • M XYZ→RGB •    Λ L D65 (λ) • S(λ) • x(λ)dλ Λ L D65 (λ) • S(λ) • ȳ(λ)dλ Λ L D65 (λ) • S(λ) • z(λ)dλ    , (8)
where

• Y D65 = Λ L D65 (λ) • ȳ(λ)dλ , • S(λ) is the reflective spectrum.    R G B    = 1 Y D65 • M XYZ→RGB •    Λo Λ i L D65 (λ i ) • S(λ i , λ o ) • x(λ o )dλ i dλ o Λo Λ i L D65 (λ i ) • S(λ i , λ o ) • ȳ(λ o )dλ i dλ o Λo Λ i L D65 (λ i ) • S(λ i , λ o ) • z(λ o )dλ i dλ o    . (9)

Using a Different Conversion Pipeline

One can use one's own conversion pipeline for a specific application. In this case, additional information can be provided as OpenEXR header attributes: the XYZ sensitivity curves, the illuminant spectrum for reflective images, and an additional exposure compensation value (see Table 6). This additional information is especially important for bi-spectral images where the illuminant greatly influences the RGB conversion even for two illuminants with the same white point. The exposition compensation is defined as follows:

RGB = RGB • 2 EV , ( 10
)
where EV is the exposure compensation value (unitless). Optional metadata for the spectral-to-tristimulus conversion pipeline. The spectrum format follows the same convention as the one described in Section 4.2.5.

Description

Specification of the full conversion pipeline, from spectral to linear RGB values, is necessary to validate one's specific conversion pipeline. With this information on hand, a producer and a consumer of a spectral image with RGB channels can validate the consistency of their spectral conversion pipelines.

Conclusion

We introduce a versatile, unambiguous, and flexible layout for storing spectral data within the existing framework of the OpenEXR image format. This layout supports emissive and reflective spectral images, and also handles polarization and bi-spectral reflective spectra. In addition to pixel data, we describe how to retain metadata from capture devices: the camera spectral response, lens transmission, and filter properties. It also allows the storage of all spectral and capture information in a single file, and common graphic software can display the spectral image in its tristimulus form if the RGB channels are also provided. The format also intrinsically supports various compression schemes on the individual channels, by virtue of being part of the OpenEXR framework. Finally, the format is modular and allows additional metadata and layers for specific applications.

0

  dB: 10 • log(φ(λ i , λ o )) (a) A full re-radiation matrix. efficiency λ o [nm] (b) Re-emission for λ i = 550 nm.

Figure 2 .

 2 Figure 2. (a) A full re-radiation matrix like this example completely describes bi-spectral reflectance for the given input (λ i ) and output (λ o ) wavelength ranges. For a given illumination wavelength λ i , the light re-emitted at various wavelengths is shown as a row in the 2D diagram. (b) A plot of the row highlighted in purple in the full matrix at λ i = 550 nm, including the diagonal part corresponding to the non-bi-spectral reflectance.

  layer = "((S[0-3])|T)" value = "(\d * ,?\d * ([eE][-+]?\d+)?

  Varying sampling rate over the spectrum.

Figure 4 .

 4 Figure 4. Default spectral bin filtering behavior when no filter is specified in the header.

  620.5nm:0;650nm:1;679.5nm:0; 620.5 650 679.5 1 [nm] Vol. 10, No. 3, 2021 http://jcgt.org (a) GIMP. (b) Spectral Viewer.

Figure 5 .

 5 Figure 5. The same spectral OpenEXR image opened in two different programs: first, a standard graphic software, GIMP [GIMP 1996], which is only capable of handling RGB images, and second, our specialized software for spectral images, Spectral Viewer [Dufay et al. 2018]. (a) Gimp cannot access the spectral information but still can render the RGB version thanks to the additional R, G, and B channels provided in the spectral OpenEXR image. (b) Our specialized software can read the spectral information stored in each pixel.

2 .

 2 Project the spectrum to CIE-XYZ 2 • 1931 color space. 3. Normalize the components with CIE D65 Y component (CIE 1931 Y).

4 .

 4 Convert XYZ to sRGB D65 linear.

Table 1 .

 1 Description of channel naming convention.

	4.1.1. Layers
	Emissive images. Layers that store individual Stokes components of a polarized
	image are named S0, S1, S2, and S3. If there is no polarization information, only
	the S0 layer is present.

Reflective images. Layers that store reflective or attenuation values are named T.

4.1.2. Channels

Each channel is named according to the wavelength or frequency bracket it represents. Because the OpenEXR format reserves dots for separating different layers, we need to use a decimal comma notation for wavelength or frequency values. Wave-

Table 2 .

 2 Mandatory version information of the layout in the OpenEXR header.

Table 4 .

 4 Mandatory metadata in the OpenEXR header for emissive polarized spectral images.

	Vol. 10, No. 3, 2021
	http://jcgt.org

Table 6 .

 6

		Attribute Name	Type	Default
	Illuminant	illuminant	Spectrum CIE D65
	X color matching function X	Spectrum CIE X 1931 2 •
	Y color matching function Y	Spectrum CIE Y 1931 2 •
	Z color matching function Z	Spectrum CIE Z 1931 2 •
	Exposure compensation	EV	Float	0

The corresponding photometric unit is nits, which is used when working with such quantities in a tristimulus pipeline.

At least mostly so: strongly polarized light is weakly visible to some observers, but the effect is not thought to significantly impact how we perceive scenes in general.

Vol. 10, No. 3, 2021 http://jcgt.org

Value must be set to 1.0 for the layout described in this article.

Only for polarized images.

Only for emissive images. Value must be "W," "W.mˆ-2," "W.srˆ-1," or "W.mˆ-2.srˆ-1."

Comma separated integer and decimal parts.

Acknowledgements

We want to thank the reviewers for their meaningful advice, among these, contacting the OpenEXR developers. We also want to thank the OpenEXR developers and mailing list members for their valuable feedback. This work was funded in part by French National Research Agency grant ANR-17-CE23-0017 "VIDA." We also acknowledge funding by Czech Science Foundation grant 19-07626S, EU H2020-MSCA-ITN-2020 grant number 956585 (PRIME), and Charles University grant number SVV-260588.

An OpenEXR Layout for Spectral Images Vol. 10, No. 3,2021 http://jcgt.org (a) When looking toward the light source, observe left circular polarized light. We refer to this as the "left"-handed convention.

(b) Looking along the light ray, we observe right circular polarized light in the same scenario as before. This is what we call the "right"-handed convention.

Capture Information

We use optional custom attributes to store properties of the acquisition device (see Table 5):

• the lens system transmission spectrum,

• the camera sensitivity curve,

• a response curve for each S0 channel.

In the case that the lens transmission spectrum or the camera sensitivity spectrum metadata field is not present, we assume that this optical element is perfect. When no filter is specified in the header for a spectral bin, the filter is assumed to be a gate function. The filter boundaries [B n (λ -), B n (λ + )] for a given bin wavelength B n are

Description

Attribute Name Type Lens system transmission spectrum lensTransmission Spectrum Camera sensitivity spectrum cameraResponse Spectrum Filter/channel bandwidth Same as S0 channel names Spectrum

Index of Supplemental Materials

Sample Code

We provide C++ sample code for reading and writing spectral and bi-spectral OpenEXR files under BSD License 2.0 terms.

• Spectral: The spectral classes ([EXR]?SpectralImage. * ) are a simpler version for one interested only in the spectral part of the OpenEXR image. This reads as well bi-spectral files but ignores the re-radiation information.

• Bi-spectral: The bi-spectral classes ( * BiSpectralImage. * ) support both spectral and bi-spectral OpenEXR images.

Sample Images

We provide several sample images.

Emissive Images

• D65.exr: The spectral power distribution of the CIE D65 illuminant. A 1 × 1 pixel image giving an example of the application for storing spectra. • BiSpectral.exr: A bi-spectral texture made of two different fluorescent reflectances.


Gilles Villard
email: gilles.villard@ens-lyon.fr

Kaltofen's division-free determinant algorithm differentiated for matrix adjoint computation

Keywords: matrix determinant, matrix adjoint, matrix inverse, characteristic polynomial, exact algorithm, division-free complexity, Wiedemann algorithm, automatic differentiation

come

Introduction

Kaltofen has proposed in [START_REF] Kaltofen | On computing determinants without divisions[END_REF] a new approach for computing matrix determinants. This approach has brought breakthrough ideas for improving the complexity estimate for the problem of computing the determinant without divisions over an abstract ring (see [START_REF] Kaltofen | On computing determinants without divisions[END_REF][START_REF] Kaltofen | On the complexity of computing determinants[END_REF]). With these foundations, the algorithm of [START_REF] Kaltofen | On the complexity of computing determinants[END_REF] computes the determinant in O(n 2.7 ) additions, subtractions, and multiplications. The same ideas also lead to the currently best known bit complexity estimate of [START_REF] Kaltofen | On the complexity of computing determinants[END_REF] for the problem of computing the characteristic polynomial.

We consider the straigth-line programs of [START_REF] Kaltofen | On computing determinants without divisions[END_REF] for computing the determinant over abstract fields or rings (with or without divisions). Using the reverse mode of automatic differentiation (see [START_REF] Linnainmaa | The representation of the cumulative rounding error of an algorithm as a Taylor expansion of the local rounding errors[END_REF][START_REF] Linnainmaa | Taylor expansion of the accumulated rounding errors[END_REF], and [START_REF] Ostrowski | Über die Berechnung von Ableitungen (in German)[END_REF])), a straight-line program for computing the determinant of a matrix A can be (automatically) transformed into a program for computing the adjoint matrix A * of A. This principle, stated by Baur and Strassen (1983, Cor. 5), is also applied by Kaltofen (1992, Sec. 1.2) for computing A * . Since the adjoint program is derived by an automatic process, few is known about the way it computes the adjoint. The only available information seems to be the determinant program itself, and the knowledge we have on the differentiation process. Neither the adjoint program can be described, or implemented, without resorting to an automatic differentiation tool.

In this paper, by studying the differentiation of Kaltofen's determinant algorithm step by step, we produce an "explicit" adjoint algorithm. The determinant algorithm, that we first recall in Section 2 over an abstract field K, uses a Krylov subspace construction, hence mainly reduces to vector times matrix, and matrix times matrix products. Another operation involved is computing the minimum polynomial of a linearly generated sequence. We apply the program differentiation mechanism, reviewed in Section 3, to the different steps of the determinant program in Section 4. This leads us to the description of a corresponding new adjoint program over a field, in Section 5. The algorithm we obtain somehow calls to mind the matrix factorization of Eberly (1997, (3.4)). We note that our objectives are similar to Eberly's ones, whose question was to give an explicit inversion algorithm from the parallel determinant algorithm of [START_REF] Kaltofen | Processor efficient parallel solution of linear systems over an abstract field[END_REF].

Our motivation for studying the differentiation and resulting adjoint algorithm, is the importance of the determinant approach of [START_REF] Kaltofen | On computing determinants without divisions[END_REF], and [START_REF] Kaltofen | On the complexity of computing determinants[END_REF], for various complexity estimates. Recent advances around the determinant of polynomial or integer matrices (see [START_REF] Eberly | Computing the determinant and Smith form of an integer matrix[END_REF]; [START_REF] Kaltofen | On the complexity of computing determinants[END_REF]; [START_REF] Storjohann | The shifted number system for fast linear algebra on integer matrices[END_REF][START_REF] Storjohann | The shifted number system for fast linear algebra on integer matrices[END_REF]), and matrix inversion (see [START_REF] Jeannerod | Asymptotically fast polynomial matrix algorithms for multivariable systems[END_REF][START_REF] Jeannerod | Asymptotically fast polynomial matrix algorithms for multivariable systems[END_REF][START_REF] Storjohann | On the complexity of inverting integer and polynomial matrices[END_REF]) also justify the study of the general adjoint problem.

For computing the determinant without divisions over a ring R, Kaltofen applies the avoidance of divisions of [START_REF] Strassen | Vermeidung von Divisionen[END_REF] to his determinant algorithm over a field. We apply the same strategy for the adjoint. From the algorithm of Section 5 over a field, we deduce an adjoint algorithm over an arbitrary ring R in Section 6. The avoidance of divisions involves computations with truncated power series. A crucial point in Kaltofen's approach is a "baby steps/giant steps" scheme for reducing the corresponding power series arithmetic cost. However, since we use the reverse mode of differentiation, the flow of computation is modified, and the benefit of the baby steps/giant steps is partly lost for the adjoint. This asks us to introduce an early, and lazy polynomial evaluation strategy for not increasing the complexity estimate.

The division-free determinant algorithm of [START_REF] Kaltofen | On computing determinants without divisions[END_REF] uses O˜(n 3.5 ) operations in R. The adjoint algorithm we propose has essentially the same cost. Our study may be seen as a first step for the differentiation of the more efficient algorithm of [START_REF] Kaltofen | On the complexity of computing determinants[END_REF]. The latter would require, in particular, to consider asymptotically fast matrix multiplication algorithms that are not discussed in what follows.

Especially in our matrix context, we note that interpreting programs obtained by automatic differentiation, may have connections with the interpretation of programs derived using the transposition principle. We refer for instance to the discussion of Kaltofen (2000, Sec. 6).

Cost functions. We let M(n) be such that two univariate polynomials of degree n over an arbitrary ring R can be multiplied using M(n) operations in R. The algorithm of [START_REF] Cantor | On fast multiplication of polynomials over arbitrary algebras[END_REF] allows

M(n) = O(n log n log log n). The function O(M(n))
also measures the cost of truncated power series arithmetic over R. For bounding the cost of polynomial gcd-type computations over a commutative field K we define the function G. Let G(n) be such that the extended gcd problem (see (von zur Gathen and Gerhard, 1999, Chap. 11)) can be solved with G(n) operations in K for polynomials of degree 2n in K[x]. The recursive Knuth/Schönhage half-Gcd algorithm (see [START_REF] Knuth | The analysis of algorithms[END_REF][START_REF] Schönhage | Schnelle Berechnung von Kettenbruchenwicklungen[END_REF][START_REF] Moenck | Fast computation of Gcds[END_REF]

) allows G(n) = O(M(n) log n).
The minimum polynomial of degree n, of a linearly generated sequence given by its first 2n terms, can be computed in (von zur Gathen and Gerhard, 1999, Algorithm 12.9)). We will often use the notation O˜that indicates missing factors of the form α(log n) β , for two positive real numbers α and β.

G(n) + O(n) operations (see

Kaltofen's determinant algorithm over a field

Kaltofen's determinant algorithm extends the Krylov-based method of [START_REF] Wiedemann | Solving sparse linear equations over finite fields[END_REF]. The latter approach is successful in various situations. We refer especially to the algorithms of [START_REF] Kaltofen | Processor efficient parallel solution of linear systems over an abstract field[END_REF] and [START_REF] Kaltofen | On Wiedemann's method of solving sparse linear systems[END_REF] around exact linear system solution that has served as basis for subsequent works. We may also point out the various questions investigated by [START_REF] Chen | Efficient matrix preconditioners for black box linear algebra[END_REF], and references therein.

Let K be a commutative field. We consider A ∈ K n×n , u ∈ K 1×n , and v ∈ K n×1 . We introduce the Hankel matrix H = uA i+j-2 v 1≤i,j≤n ∈ K n×n , and let h k = uA k v for 0 ≤ k ≤ 2n -1. We also assume that H is non-singular:

det H = det         uv uAv . . . uA n-1 v uAv uA 2 v . . . uA n v . . . . . . . . . . . . uA n-1 v . . . . . . uA 2n-2 v         = 0. ( 1
)
In the applications, (1) is ensured either by construction of A, u, and v, as in [START_REF] Kaltofen | On computing determinants without divisions[END_REF][START_REF] Kaltofen | On the complexity of computing determinants[END_REF], or by randomization (see the above cited references around Wiedemann's approach, and [START_REF] Kaltofen | On computing determinants without divisions[END_REF][START_REF] Kaltofen | On the complexity of computing determinants[END_REF]).

One of the key ideas of [START_REF] Kaltofen | On computing determinants without divisions[END_REF] for reducing the division-free complexity estimate for computing the determinant, is to introduce a "baby steps/giant steps" behaviour in the Krylov subspace construction. With baby steps/giant steps parameters r = ⌈2n/s⌉ and s = ⌈ √ n⌉ (rs ≥ 2n) we consider the following algorithm.

Algorithm Det

Input: A ∈ K n×n , u ∈ K 1×n , v ∈ K n×1 step i. v 0 := v; For i = 1, . . . , r -1 do v i := Av i-1 step ii. B := A r step iii. u 0 := u; For j = 1, . . . , s -1 do u j := u j-1 B step iv. For i = 0, 1, . . . , r -1 do For j = 0, 1, . . . , s -1 do h i+jr := u j v i step v. f := the minimum polynomial of {h k } 0≤k≤2n-1 Output: det A := (-1) n f (0).
We ommit the proof of next theorem that establishes the correctness and the cost of Algorithm Det, and refer to [START_REF] Kaltofen | On computing determinants without divisions[END_REF]. We may simply note that the sequence {h k } 0≤k≤2n-1 is linearly generated. In addition, if (1) is true, then the minimum polynomial f of {h k } 0≤k≤2n-1 , the minimum polynomial of A, and the characteristic (monic) polynomial of A coincide. Hence (-1) n f (0) is equal to the determinant of A. Via an algorithm that can multiply two matrices of K n×n in O(n ω ) we have:

Theorem 1. If A ∈ K n×n , u ∈ K 1×n , and v ∈ K n×1 satisfy (1), then Algorithm Det computes the determinant of A in O(n ω log n) operations in K.
For the matrix product we may set ω = 3, or ω = 2.376 using the algorithm of [START_REF] Coppersmith | Matrix multiplication via arithmetic progressions[END_REF]. In the rest of the paper we work with a cubic matrix multiplication algorithm. Our study has to be generalized if fast matrix multiplication is introduced.

Backward automatic differentiation

The determinant of A ∈ K n×n is a polynomial ∆ in K[a 1,1 , . . . , a i,j , . . . , a n,n ] of the entries of A. We denote the adjoint matrix by A * such that AA * = A * A = (det A)I. As noticed by [START_REF] Baur | The complexity of partial derivatives[END_REF], the entries of A * satisfy

a * j,i = ∂∆ ∂a i,j , 1 ≤ i, j ≤ n. (2)
The reverse mode of automatic differentiation allows to transform a program which computes ∆ into a program which computes all the partial derivatives in (2). Among the rich literature about the reverse mode of automatic differentiation we may refer to the seminal works of [START_REF] Linnainmaa | The representation of the cumulative rounding error of an algorithm as a Taylor expansion of the local rounding errors[END_REF][START_REF] Linnainmaa | Taylor expansion of the accumulated rounding errors[END_REF] and [START_REF] Ostrowski | Über die Berechnung von Ableitungen (in German)[END_REF]. For deriving the adjoint program from the determinant program we follow the lines of [START_REF] Baur | The complexity of partial derivatives[END_REF] and [START_REF] Morgenstern | How to compute fast a function and all its derivatives, a variation on the theorem of Baur-Strassen[END_REF].

Algorithm Det is a straight-line program over K. For a comprehensive study of straight-line programs for instance see (Bürgisser et al., 1997, Chapter 4). We assume that the entries of A are stored initially in n 2 variables δ i , -n 2 < i ≤ 0. Then we assume that the algorithm is a sequence of arithmetic operations in K, or assignments to constants of K. Let L be the number of such operations. We assume that the result of each instruction is stored in a new variable δ i , hence the algorithm is seen as a sequence of instructions

δ i := δ j op δ k , op ∈ {+, -, ×, ÷}, -n 2 < j, k < i, (3)
or

δ i := c, c ∈ K, (4) for 1 ≤ i ≤ L.
Note that a binary arithmetic operation (3) where one of the operands is a constant of K can be implemented with the aid of (4). For any 0 ≤ i ≤ L, the determinant maybe be seen as a rational function

∆ i of δ -n 2 +1 , . . . , δ i , such that ∆ 0 (δ -n 2 +1 , . . . , δ 0 ) = ∆(a 1,1 , . . . , a n,n ), (5)
and such that the last instruction gives the result:

det A = δ L = ∆ L (δ -n 2 +1 , . . . , δ L ). ( 6
)
The reverse mode of automatic differentiation computes the derivatives (2) in a backward recursive way, from the derivatives of ( 6) to those of ( 5). Using ( 6) we start the recursion with

∂∆ L ∂δ L = 1, ∂∆ L ∂δ l = 0, -n 2 < l ≤ L -1.
Then, writing

∆ i-1 (δ -n 2 +1 , . . . , δ i-1 ) = ∆ i (δ -n 2 +1 , . . . , δ i ) = ∆ i (δ -n 2 +1 , . . . , g(δ j , δ k )), (7)
where g is given by ( 3) or (4), we have

∂∆ i-1 ∂δ l = ∂∆ i ∂δ l + ∂∆ i ∂δ i ∂g ∂δ l , -n 2 < l ≤ i -1, (8)
for 1 ≤ i ≤ L. Depending on g several cases may be examined. For instance, for an addition δ i := g(δ k , δ j ) = δ k + δ j , (8) becomes

∂∆ i-1 ∂δ k = ∂∆ i ∂δ k + ∂∆ i ∂δ i , ∂∆ i-1 ∂δ j = ∂∆ i ∂δ j + ∂∆ i ∂δ i , (9)
with the other derivatives (l = k or j) remaining unchanged. In the case of a multiplication (8) gives that the only derivatives that are modified are

δ i := g(δ k , δ j ) = δ k × δ j ,
∂∆ i-1 ∂δ k = ∂∆ i ∂δ k + ∂∆ i ∂δ i δ j , ∂∆ i-1 ∂δ j = ∂∆ i ∂δ j + ∂∆ i ∂δ i δ k . ( 10
)
We see for instance in (10), where δ j is used for updating the derivative with respect to δ k , that the recursion uses intermediary results of the determinant algorithm. For the adjoint algorithm, we will assume that the determinant algorithm has been executed once, and that the δ i 's are stored in n 2 + L memory locations.

Recursion (8) gives a practical mean, and a program, for computing the N = n 2 derivatives of ∆ with respect to the a i,j 's. For any rational function Q in N variables δ -N +1 , . . . , δ 0 the corresponding general statement is: Theorem 2. [START_REF] Baur | The complexity of partial derivatives[END_REF]] Let P be a straight-line program computing Q in L operations in K. One can derive an algorithm ∂P that computes Q and the N partial derivatives ∂Q/∂δ l in less than 5L operations in K.

Combining Theorem 2 with Theorem 1 gives the construction of an algorithm ∂Det for computing the adjoint matrix A * (see (Baur and Strassen, 1983, Corollary 5)). The algorithm can be generated automatically via an automatic differentiation tool1 . However, it seems unclear how it could be programmed directly, and, to our knowledge, it has no interpretation of its own.

Differentiating the determinant algorithm over a field

We apply the backward recursion (8) to Algorithm Det of Section 2 for deriving the algorithm ∂Det. We assume that A is non-singular, hence A * is non-trivial. By construction, the flow of computation for the adjoint is reversed compared to the flow of Algorithm Det, therefore we start with the differentiation of step v.

Differentiation of the minimum polynomial constant term computation

At step v, Algorithm Det computes the minimum polynomial f of the linearly generated sequence {h k } 0≤k≤2n-1 . Let λ be the first instruction index at which all the h k 's are known. We apply the recursion until step λ, globally, we mean that we compute the derivatives of ∆ λ . After the instruction λ, the determinant is viewed as a function ∆ v of the h k 's only. Following (7) we have

det(A) = ∆ λ (δ -n 2 +1 , . . . , δ λ ) = ∆ v (h 1 , . . . , h 2n-1 ).
Hence we may focus on the derivatives ∂∆ v /∂h k , 0 ≤ k ≤ 2n -1, the remaining ones are zero.

Using assumption (1) we know that the minimum polynomial f of {h k } 0≤k≤2n-1 has degree n, and if

f (x) = f 0 + f 1 x + . . . + f n-1 x n-1 + x n , then f satisfies H         f 0 f 1 . . . f n-1         =         h 0 h 1 . . . h n-1 h 1 h 2 . . . h n . . . . . . . . . . . . h n-1 . . . . . . h 2n-2                 f 0 f 1 . . . f n-1         = -         h n h n+1 . . . h 2n-1         (11)
see, e.g., [START_REF] Kaltofen | On computing determinants without divisions[END_REF], or (von zur Gathen and Gerhard, 1999, Algorithm 12.9) together with [START_REF] Brent | Fast solution of Toeplitz systems of equations and computation of Padé approximations[END_REF]. Applying Cramer's rule we see that

f 0 = (-1) n det         h 1 h 2 . . . h n h 2 h 3 . . . h n+1 . . . . . . . . . . . . h n . . . . . . h 2n-1         / det H, hence, defining H A = uA i+j-1 v 1≤i,j≤n = (h i+j-1 ) 1≤i,j≤n ∈ K n×n , we obtain ∆ v = det H A det H . ( 12
)
Let Ku and K v be the Krylov matrices

Ku = [u T , A T u T , . . . , (A T ) n-1 u T ] T ∈ K n×n , (13)
and

K v = [v, Av, . . . , A n-1 v] ∈ K n×n . ( 14
) Since H = Ku K v , assumption (1) implies that both Ku and K v are non-singular. Hence, using that A is non-singular, we note that H A = Ku AK v also is non-singular.

For differentiating ( 12), let us first specialize (2) to Hankel matrices. We denote by (∂∆/∂a i,j )(H) the substitution of the a i,j 's for the entries of H in ∂∆/∂a i,j , for 1 ≤ i, j ≤ n. From (2) we have

h * j,i = ∂∆ ∂a i,j (H), 1 ≤ i, j ≤ n.
Since the entries of H are constant along the anti-diagonals, we deduce that

∂ det H ∂h k = i+j-2=k ∂∆ ∂a i,j (H) = i+j-2=k h * j,i = i+j-2=k h * i,j , 0 ≤ k ≤ 2n -2.
In other words, we may write

∂ det H ∂h k = σ k (H * ), 0 ≤ k ≤ 2n -1, (15)
where, for a matrix M = (m ij ), we define

σ k (M ) = 0 + i+j-2=k m ij , 1 ≤ i, j ≤ n.
The function σ k (M ) is the sum of the entries in the anti-diagonal of M starting with

m 1,k+1 if 0 ≤ k ≤ n -1, and m k-n+2,n if n ≤ k ≤ 2n -2.
Shifting the entries of H for obtaining H A we also have

∂ det H A ∂h k = σ k-1 (H * A ), 0 ≤ k ≤ 2n -1. ( 16
)
Now, differentiating (12), together with ( 15) and ( 16), leads to

∂∆ v ∂h k = (∂ det H A /∂h k ) det H - (∂ det H/∂h k ) det H det H A det H = (∂ det H A /∂h k ) det H A det H A det H -σ k (H -1 )∆ v
and, consequently, to

∂∆ v ∂h k = σ k-1 (H -1 A ) -σ k (H -1 ) ∆ v , 0 ≤ k ≤ 2n -1. ( 17
)
With ( 17) we identify the problem solved by the first step of the ∂Det algorithm, and provide first informations for interpreting or implementing the adjoint program. Various algorithms may be used for computing the minimum polynomial (for instance see (von zur Gathen and Gerhard, 1999, Algorithm 12.9)), that will lead to corresponding algorithms for computing the left sides in (17). However, we will not discuss these aspects, since the associated costs are not dominant in the overall complexity.

We have recalled, in the introduction, that the minimum polynomial f (its constant term f (0)) can be computed from the h k 's in G(n) + O(n) operations in K. Hence Theorem 2 gives an algorithm for computing the derivatives using 5G(n) + O(n) operations.

Alternatively, in the Appendix we propose a direct approach that takes advantage of (17). Proposition 4 shows that if f , H, and H A are given, then the ∂∆ v /∂h k 's can be computed in G(n) + O(M(n)) operations in K.

Differentiation of the dot products

For differentiating step iv, ∆ is seen as a function ∆ iv of the u j 's and v i 's. The entries of u j are used for computing the r scalars h jr , h 1+jr , . . . , h (r-1)+jr for 0 ≤ j ≤ s -1. The entries of v i are involved in the computation of the s scalars h i , h i+r , . . . , h i+(s-1)r for 0 ≤ i ≤ r -1.

In ( 8), the new derivative ∂∆ i-1 /∂δ l is obtained by adding the current instruction contribution to the previously computed derivative ∂∆ i /∂δ l . Since all the h i+jr 's are computed independently according to

h i+jr = n l=0 (u j ) l (v i ) l ,
it follows that the derivative of ∆ iv with respect to an entry (u j ) l or (v i ) l is obtained by summing up the contributions of the multiplications (u j ) l (v i ) l . We obtain

∂∆ iv ∂(u j ) l = r-1 i=0 ∂∆ v ∂h i+jr (v i ) l , 0 ≤ j ≤ s -1, 1 ≤ l ≤ n, (18)
and

∂∆ iv ∂(v i ) l = s-1 i=0 ∂∆ v ∂h i+jr (u j ) l , 0 ≤ i ≤ r -1, 1 ≤ l ≤ n. ( 19
)
By abuse of notations (of the sign ∂), we let ∂u j be the n × 1 vector, respectively ∂v i be the 1 × n vector, whose entries are the derivatives of ∆ iv with respect to the entries of u j , respectively v i . Note that because of the index transposition in (2), it is convenient, here and in the following, to take the transpose form (column versus row) for the derivative vectors. Defining also

∂H = ∂∆ v ∂h i+jr 0≤i≤r-1, 0≤j≤s-1 ∈ K r×s ,
we deduce, from ( 18) and ( 19), that

[∂u 0 , ∂u 1 , . . . , ∂u s-1 ] = [v 0 , v 1 , . . . , v r-1 ] ∂H ∈ K n×s . ( 20
)
and

        ∂v 0 ∂v 1 . . . ∂v r-1         = ∂H         u 0 u 1 . . . u s-1         ∈ K r×n . (21)
Identities ( 20) and ( 21) give the second step of the adjoint algorithm. In Algorithm Det, step iv costs essentially 2rsn additions and multiplications in K. Here we have essentially 4rsn additions and multiplications using basic loops (as in step iv) for calculating the matrix products, we mean without an asymptotically fast matrix multiplication algorithm.

Differentiation of the matrix times vector and matrix products

The recursive process for differentiating step iii to step i may be written in terms of the differentiation of the basic operation (or its transposed operation)

q := p • M ∈ K 1×n , (22)
where p and q are row vectors of dimension n, and M is an n × n matrix. We assume at this point (by construction of the recursion) that column vectors ∂p and ∂q of derivatives of the determinant with respect to the entries of p and q, are available. For instance, for differentiating step iii, we will consider the ∂u j 's. We also assume that an n × n matrix ∂M , whose transpose gives the derivatives with respect to the m ij 's, has been computed. Initially, for step iii, we will take ∂B = 0. Following the lines of previous section for obtaining ( 20) and ( 21), we see that differentiating ( 22) amounts to updating ∂p and ∂M according to

   ∂p := ∂p + M • ∂q ∈ K n , ∂M := ∂M + ∂q • p ∈ K n×n . (23)
Starting from the values of the ∂u j 's computed with (20), and from ∂B = 0, for the differentiation of step iii, (23) gives

   ∂u j-1 := ∂u j-1 + B • ∂u j , ∂B := ∂B + ∂u j • u j-1 , j = s -1, . . . , 1. (24)
For step ii, we mean B := A r , we show that the backward recursion leads to

∂A := r k=1 A r-k • ∂B • A k-1 . ( 25
)
Here, the notation ∂A stands for the n × n matrix whose transpose gives the derivatives ∂∆ ii /∂a i,j . We may show (25) by induction on r. For r = 1, ∂A = ∂B is true. If ( 25) is true for r -1, then let C = A r-1 and B = CA. Using (23), and overloading the notation ∂A, we have

   ∂C = A • ∂B ∈ K n×n , ∂A = ∂B • C ∈ K n×n .
Hence, using (25) for r -1, we establish that

∂A = ∂A + r-1 k=1 A r-k-1 • ∂C • A k-1 , = ∂B • C + r-1 k=1 A r-k-1 • (A • ∂B) • A k-1 = ∂B • A r-1 + r-1 k=1 A r-k • ∂B • A k-1 = r k=1 A r-k • ∂B • A k-1
. Any specific approach for computing A r will lead to an associated program for computing ∂A. Let us look, in particular, at the case where step ii of Algorithm Det is implemented by repeated squaring, in essentially log 2 r matrix products. Consider the recursion

A 0 := A For k = 1, . . . , log 2 r do A 2 k := A 2 k-1 • A 2 k-1 B := A r
that computes B := A r . The associated program for computing the derivatives is

∂A r := ∂B For k = log 2 r, . . . , 1 do ∂A 2 k-1 := A 2 k-1 • ∂A 2 k + ∂A 2 k • A 2 k-1 ∂A := ∂A 0 , (26)
and costs essentially 2 log 2 r matrix products.

From the values of the ∂v i 's computed with (21), we finally differentiate step i, and update ∂A according to

   ∂v i-1 := ∂v i-1 + ∂v i • A, ∂A := ∂A + v i-1 • ∂v i , i = r -1, . . . , 1. (27)
Now, ∂A is the n × n matrix whose transpose gives the derivatives ∂∆ i /∂a i,j = ∂∆/∂a i,j , hence from (2) we know that A * = ∂A.

step iii and step i both cost essentially r (≈ s) matrix times vector products. From ( 24) and ( 27) the differentiated steps both require r matrix times vector products, and 2rn 2 + O(rn) additional operations in K.

The adjoint algorithm over a field

We call Adjoint the algorithm obtained from the successive differentiations of Section 4. Algorithm Adjoint is detailed below. We keep the notations of previous sections. We use in addition U ∈ K s×n and V ∈ K n×r (resp. ∂U ∈ K n×s and ∂V ∈ K r×n ) for the right sides (resp. the left sides) of ( 20) and ( 21).

The cost of Adjoint is dominated by step iv * , which is the differentiation of the matrix power computation. As we have seen with ( 26), the number of operation is essentially twice as much as for Algorithm Det. The code we give allows an easy implementation.

We note that if the product by det A is avoided in step i * , then the algorithm computes the matrix inverse A -1 . We may put this into perspective with the algorithm given by [START_REF] Eberly | Processor-efficient parallel matrix inversion over abstract fields: two extensions[END_REF]. With Ku and K v the Krylov matrices of ( 13) and ( 14), Eberly has proposed a processor-efficient inversion algorithm based on

A -1 = K v H -1 A Ku . ( 28
)
To see whether a baby steps/giant steps version of ( 28) would lead to an algorithm similar to Adjoint deserves further investigations.

Algorithm Adjoint (∂Det)

Input: A ∈ K n×n non-
:= σ k-1 (H -1 A ) -σ k (H -1 ) det A, 0 ≤ k ≤ 2n -1 step ii * .
/* Requires the u j 's and v i 's, see ( 20) and ( 21

) */ ∂U := V • ∂H ∂V := ∂H • U step iii * . /* Requires B = A r , see (24) */ For j = s -1, . . . , 1 do ∂u j-1 := ∂u j-1 + B • ∂u j ∂B := ∂B + ∂u j • u j-1
step iv * . /* Requires the powers of A, see ( 25) or ( 26) */

A * := r k=1 A r-k • ∂B • A k-1 step v * . /* See (27) */ For i = r -1, . . . , 1 do ∂v i-1 := ∂v i-1 + ∂v i • A A * := A * + v i-1 • ∂v i
Output: The adjoint matrix A * ∈ K n×n .

Application to computing the adjoint without divisions

Now let A be an n × n matrix over an abstract ring R. Kaltofen's algorithm for computing the determinant of A without divisions applies Algorithm Det on a well chosen univariate polynomial matrix Z(z) = C + z(A -C) where C ∈ Z n×n , with a dedicated choice of projections u = ϕ ∈ Z 1×n and v = ψ ∈ Z n×1 . The algorithm uses Strassen's avoidance of divisions (see [START_REF] Strassen | Vermeidung von Divisionen[END_REF][START_REF] Kaltofen | On computing determinants without divisions[END_REF]). Since the determinant of Z is a polynomial of degree n in z, the arithmetic operations over K in Det may be replaced by operations on power series in R

[[z]] modulo z n+1 . Once the determinant of Z(z) is computed, the evaluation (det Z)(1) = det(C + 1 × (A -C))
gives the determinant of A. The choice of C, ϕ and ψ is such that, whenever a division by a truncated power series is performed the constant coefficients are ±1. Therefore the algorithm necessitates no divisions. Note that, by construction of Z(z), the constant terms of the power series involved when Det is called with inputs Z(z), ϕ and ψ, are the intermediary values computed by Det with inputs C, ϕ and ψ.

The cost for computing the determinant of A without divisions is then deduced as follows. In step i and step ii of Algorithm Det applied to Z(z), the vector and matrix entries are polynomials of degree O( √ n). The cost of step ii dominates, and is

O(n 3 M( √ n) log n) = O˜(n 3 √ n) operations in R. step iii, iv, and v cost O(n 2 √ n) op- erations on power series modulo z n+1 , that is O(n 2 M(n) √ n) operations in R. Hence det Z(z) is computed in O˜(n 3 √ n) operations in R,
and det A is obtained with the same cost bound.

An main property of Kaltofen's approach (which also holds for the improved blocked version of [START_REF] Kaltofen | On the complexity of computing determinants[END_REF]), is that the scalar value det A is obtained via the computation of the polynomial value det Z(z). This property seems to be lost with the adjoint computation. We are going to see how Algorithm Adjoint applied to Z(z) allows to compute A * ∈ R n×n in time O˜(n 3 √ n) operations in R, but does not seem to allow the computation of Z * (z) ∈ R[z] n×n with the same complexity estimate. Indeed, a key point in Kaltofen's approach for reducing the overall complexity estimate, is to compute with small degree polynomials (degree O( √ n)) in step i and step ii. However, since the adjoint algorithm has a reversed flow, this point does not seem to be relevant for Adjoint, where polynomials of degree n are involved from the beginning. Our approach for computing A * over R keeps the idea of running Algorithm Adjoint with input Z(z) = C + z(A -C), such that Z * (z) has degree less than n, and gives A * = Z * (1). In Section 6.1, we verify that the implementation using Proposition 4, needs no divisions. We then show in Section 6.2 how to establish the cost estimate O˜(n 3 √ n). The principle we follow is to start evaluating polynomials at z = 1 as soon as computing with the entire polynomials is prohibitive.

Division-free Hankel matrix inversion and anti-diagonal sums

In Algorithm Adjoint, divisions may only occur during the anti-diagonal sums computation. We verify here that with the matrix Z(z), and the special projections ϕ ∈ Z 1×n , ψ ∈ Z n×1 , the approach described in the Appendix for computing the anti-diagonal sums requires no divisions. Equivalently, since we use Strassen's avoidance of divisions, we verify that with the matrix C and the projections ϕ, ψ, the approach necessitates no divisions. As we are going to see, this a direct consequence of the construction of [START_REF] Kaltofen | On computing determinants without divisions[END_REF].

Here we let h k = ϕC k ψ for 0 ≤ k ≤ 2n -1, a(x) = x 2n , and b(x) = h 0 x 2n-1 + h 1 x 2n-2 + . . . + h 2n-1 . The extended Euclidean scheme with inputs a and b leads to a normal sequence, and after n -1 and n steps of the scheme, we get (see (Kaltofen, 1992, Sec. 2)):

s(x)a(x) + t(x)b(x) = c(x), with deg s = n -2, deg t = n -1, deg c = n, (29) and s
(x)a(x) + t(x)b(x) = c(x), with deg s = n -1, deg t = n, deg c = n -1. ( 30
)
The polynomial t is such that

t = ±x n + intermediate monomials + 1 = ±f, (31)
with f the minimum polynomial of {h k } 0≤k≤2n-1 . One may check, in particular, that the n equations obtained by identifying the coefficients of degree 2n -1 ≥ k ≥ n in (30) give the linear system (11), that defines f . The polynomial c also has leading coefficient ±1. By identifying the coefficients of degree 2n -1 ≥ k ≥ n in (29), we obtain:

H         t 0 t 1 . . . t n-1         =         h 0 h 1 . . . h n-1 h 1 h 2 . . . h n . . . . . . . . . . . . h n-1 . . . . . . h 2n-2                 t 0 t 1 . . . t n-1         = ±         0 0 . . . 1         . ( 32
)
Therefore t = ±g with g the polynomial needed for computing ( 44)-( 47), in addition to f . Since C, ϕ, and ψ are such that the extended Euclidean scheme necessitates no divisions (see (Kaltofen, 1992, Sec. 2)), we see that both f and g may be computed with no divisions. The only remaining division in the algorithm for Proposition 4 is at (38). From ( 31), this division is by f 0 = 1.

Lazy polynomial evaluation and division-free adjoint computation

We run Algorithm Adjoint with input Z(z) ∈ R[z] n×n , and start with operations on truncated power series modulo z n+1 . We assume that Algorithm Det has been executed, and that its intermediary results have been stored.

Using Proposition 4 and previous section,

step i * requires O(G(n)M(n)) = O˜(n 2 ) operations in R for computing ∂H(z) of degree n in R[z] r×s . step ii * , step iii * , and v * cost O(n 2 √ n) operations in K,
hence, taking into account the power series operations, this gives O(n 2 M(n)

√ n) = O˜(n 3 √ n) operations in R for the division-free version. The cost analysis of step iv * , using (26) over power series modulo z n , leads to log 2 r matrix products, hence to the time bound O˜(n 4 ), greater than the target estimate O˜(n 3 √ n). As noticed previously, step iii of Algorithm Det only involves polynomials of degree O( √ n), while the reversed program for step iv * of Algorithm Adjoint, relies on ∂B(z)

whose degree is n.

Since only Z * (1) = A * is needed, our solution, for restricting the cost to O˜(n 3 √ n), is to start evaluating at z = 1 during step iv * . However, since power series multiplications are done modulo z n , this evaluation must be lazy. The fact that matrices Z k (z), 1 ≤ k ≤ r-1, of degree at most r-1 are involved, enables the following. Let a and c be two polynomials such that deg a+deg c = r -1 in R[z], and let b be of degree n ≥ r -2 in R[z]. Considering the highest degree part of b, and evaluating the lowest degree part at z = 1, we define b

H (z) = b n z r-2 + . . . + b n-r+2 ∈ R[z] and b L = b n-r+1 + . . . + b 0 ∈ R. We then remark that a(z)b(z)c(z) mod z n+1 (1) = a(z)(b H (z)z n-r+2 + b L )c(z) mod z n+1 (1), = a(z)b H (z)c(z) mod z r-1 (1) + (a(z)b L c(z)) (1).
* := r k=1 Z r-k (z) • ∂B H (z) • Z k-1 (z) mod z r-1 (1) Z * := Z * + r k=1 Z r-k (z) • ∂B L • Z k-1 (z) (1), (34) in O˜(n 3 M(r)) = O˜(n 3 √ n) operations in R.
This leads to an intermediary value Z * ∈ R n×n before step v * . The value is updated at step v * with power series operations, and a final evaluation at z = 1 in time O˜(n 2 rM(n)) = O˜(n 3 √ n). Since only step iv * has been modified, we obtain the following result.

Theorem 3. Let A ∈ R n×n . If Algorithm Adjoint, modified according to (34), is executed with input Z(z) = C + z(A -C), power series operations modulo z n+1 , and a final evaluation at z = 1, then the matrix adjoint A * is computed in O˜(n 3 √ n) operations in R.

Concluding remarks

We have developed an explicit algorithm for computing the matrix adjoint using only ring arithmetic operations. The algorithm has complexity estimate O˜(n 3.5 ). It represents a practical alternative to previously existing solutions for the problem, that rely on automatic differentiation of a determinant algorithm. Our description of the algorithm allows direct implementations. It should help understanding how the adjoint is computed using Kaltofen's baby steps/giant steps construction. Still, a full mathematical explanation deserves to be investigated. Our work has to be generalized to the block algorithm of [START_REF] Kaltofen | On the complexity of computing determinants[END_REF] (with the use of fast matrix multiplication algorithms) whose complexity estimate is currently the best known for computing the determinant, and the adjoint without divisions.

For implementing (17), we study the computation of the anti-diagonal sums σ k of H -1 and H -1 A . We first use the formula of [START_REF] Labahn | The inverses of block Hankel and block Toeplitz matrices[END_REF] for Hankel matrices inversion. The minimum polynomial f of {h and satisfies (11). Let the last column of H -1 be given by

k } 0≤k≤2n-1 is f (x) = f 0 + f 1 x + . . . + f n-1 x n-1 + x n ,
H [g 0 , g 1 , . . . , g n-1 ] T = [0, . . . , 0, 1] T ∈ K n . (35)
Applying [START_REF] Labahn | The inverses of block Hankel and block Toeplitz matrices[END_REF], Theorem 3.1) with ( 11) and ( 35), we know that

H -1 =         f 1 . . . f n-1 1 . . . . . . . . . f n-1 . . . 0 1              g 0 . . . g n-1 . . . . . . 0 g 0      -         g 1 . . . g n-1 0 . . . . . . . . . g n-1 . . . 0 0              f 0 . . . f n-1 . . . . . . 0 f 0      . (36)
For deriving an analogous formula for H -1 A , using the notations of Section 4.1, we first recall that H = Ku K v and H A = Ku AK v . Multiplying (11) on the left by Ku A K-1 u gives

H A [f 0 , f 1 , . . . , f n-1 ] T = -[h n+1 , h n+2 , . . . , h 2n ] T . (37)
We also notice that

H A H -1 = K -1 u A T K u T ,
and, using the action of A T on the vectors u T , . . . , (A T ) n-2 u T , we check that H A H -1 is the companion matrix

H A H -1 =         0 1 0 . . . . . . 0 . . . 0 1 -f 0 -f 1 . . . -f n-1         . Hence the last column [g * 0 , g * 1 , . . . , g * n-1 ] of H -1
A is the first column of H -1 divided by -f 0 . Using (36) for determining the first column of H -1 , we get

[g * 0 , g * 1 , . . . , g * n-1 ] T = - g 0 f 0 [f 1 , . . . , f n-1 , 1] T + [g 1 , . . . , g n-1 , 0] T . (38)
Applying (Labahn et al., 1990, Theorem 3.1), now with ( 37) and ( 38), we obtain We have (43)

H -1 A =         f 1 . . .
m i,j = i+j-2 s=i-1 l s r i+j-s-2 , 1 ≤ i + j -1 ≤ n, (40)
It remains to apply ( 42) and ( 43) to the structured matrix products in ( 36) and ( 39), for computing the σ k (H -1 ) and σ k (H -1 A )'s. Together with the minimum polynomial f = f 0 +. . .+f n-1 x n-1 +x n , let g = g 0 +. . .+g n-1 x n-1 (see ( 35)), and g * = g * 0 . . .+g * n-1 x n-1 (see ( 38)). We may now combine, respectively (36) and (39), with (42), for obtaining

f ′ g -g ′ f mod x n = n-1 k=0 σ k (H -1 ) x k , (44)
and

f ′ g * -(g * ) ′ f mod x n = n-1 k=0 σ k (H -1 A ) x k . (45)
Defining also rev(f ) = 1 + f n-1 x + . . . + f 0 x n , rev(g) = g n-1 x + . . . + g 0 x n , and rev(g * ) = g * n-1 x + . . . + g * 0 x n , the combination of, respectively, (36) and (39), with (43), leads to rev(g) ′ rev(f ) -rev(f ) ′ rev(g) mod Using the approach of [START_REF] Brent | Fast solution of Toeplitz systems of equations and computation of Padé approximations[END_REF] we know that computing the last column of H -1 reduces to an extended Euclidean problem of degree 2n. Hence the polynomial g is computed in G(n) + O(n) operations. From there, g * is computed using (38). Then, applying (44)-(47) leads to the cost O(M(n)).

  step iv * , we follow the definition of b H and b L , and first compute ∂B H (z) ∈ R[z] n×n of degree r-2, and ∂B L ∈ R n×n . Applying (33), the sum r k=1 Z r-k (z)• ∂B(z) • Z k-1 (z) may then be evaluated at z = 1 by the program Modified step iv * . Z

  and (39) we see that computing σ k (H -1 ) and σ k-1 (H -1 A ), for 0 ≤ k ≤ 2n-1, reduces to computing the anti-diagonal sums for a product of triangular Hankel times triangular Toeplitz matrices. Let M = LR =

l

  s r i+j-s-2 , n ≤ i + j -1 ≤ 2n -1. (41) For 0 ≤ k ≤ 2n -2, σ k (M ) is defined by summing the m i,j 's such that i + j -2 = k. Using (40) we obtain σ k (M ) = k+1 i=1 m i,k-i+2 = k+1 i=1 k s=i-1 l s r k-s , = k s=0 (s + 1)l s r k-s , 0 ≤ k ≤ n -1 (s + n -k) l s r k-s , n -1 ≤ k ≤ 2ns-1 x s ) mod x n = n-1 k=0 σ 2n-k-2 (M ) x k .

  * ) ′ rev(f ) -rev(f ) ′ rev(g * ) mod x n = n-1 k=0 σ 2n-k-2 (H A ) x k . (47)Proposition 4. Assume that the minimum polynomial f and the Hankel matrices H and H A are given. The anti-diagonal sums σ k (H -1 ) andσ k (H -1 A ), for 0 ≤ k ≤ 2n -1, can be computed in G(n) + O(M(n)) operations in K.

  singular, and the intermediary data of Algorithm Det All the derivatives are initialized to zero step i * . /* Requires the Hankel matrices H and H A , see (17) */ ∂∆ v /∂h k

We refer for instance to http://www.autodiff.org

Acknowledgements. We thank Erich Kaltofen who has brought reference Ostrowski et al. (1971) to our attention.

Appendix: Hankel matrix inversion and anti-diagonal sums

⋆ This research was partly supported by the French National Research Agency, ANR Gecko.


Guillaume Grenet

Frédéric Alexandre

Behavior learning using emotional conditioning

We present a novel way to design a control system for a robot, using emotions as a way to produce richer internal states. We believe that using a single scalar as an evaluation of the quality of the policy and stating that the goal of the agent is to gather reward, as it is proposed by reinforcement learning, is not an appropriate granularity for creating an autonomous control system : even with a fine-tuned reward function, efficient on a specific task, it is often impractical to derive any useful knoweldge from it in order to build more flexible, neuromimetics control systems. A complete shift of paradigm is necessary for a bottom up approach : the robot is given pain and pleasure perception circuits and we examine how emotions arise and are the basis for respondant and operant conditioning. Inspired from the cerebral circuits of superior mammals responsible of this behavior, we propose an implementation in an autonomous robot using models of adaptive neural networks.

Motivations and emotions

As sensorimotor capabilities are given to a robot, the subsequent question is to wonder what makes the robot use these capabilities and act in the environment: What are its motivations? Reinforcement learning defines some states of the environment as goals/sub-goals (necessary steps to reach the main goal) and introduces a reward signal as a measure of the distance to those goals. This learning selects optimal policies from the anticipation and unfolding through time of the consequences of each possible action on the environment. Despite the impressive results achieved by reinforcement learning, it seems that it may be too naive to consider the reward signal as a proper way to describe goals : it does not promote rich internal states, as exploited by us, humans, when asked to describe our motivational strategies. This latter statement is the basis of the ongoing researches that we present in this paper: inspired by a behavioral, psychological and physiological analysis of motivation in humans, a richer model of internal states in robots could endow them with more realistic motivational characteristics. Several considerations are at the basis of the choices that have been made for the internal architecture:

1. Two different circuits related to pain and pleasure perception are the way by which we perceive states of the world (including ourselves) as goals to be avoided or looked after. Those two signals are sufficient to describe a wide range of emotions as described in [START_REF] Rolls | Precis of the brain and emotion[END_REF], for instance, fear can be seen as the anticipation of pain while frustration can be seen as the absence of an expected pleasure. While it may be argued that a positive/negative single reward signal could be enough, the use of two distinct circuits allows to deal with situations where stimuli have both positive and negative valence in a satisfactory way. It is also our purpose here to demonstrate how giving an adaptive meaning to both signals can lead to an efficient control architecture. 2. Some situations are innately "unconditioned stimuli" for which we "automatically" anticipate pain or pleasure.

New neutral situations can be associated by learning and become "conditioned stimuli". This is classical (or Pavlovian) conditioning, the rules of which have been extensively modeled [13][7]. In our approach, we propose that this learning (often associated to the limbic system) is also the basis for our discrimination capabilities (developed in the cortical temporal areas, ideally localized between the perceptual cortical regions and the limbic system): consistently with the Rescorla-Wagner rule [START_REF] Rescorla | A theory of pavlovian conditioning: Variations in the effectiveness of reinforcement and nonreinforcement[END_REF], we learn to discriminate stimuli which evoke surprise (not predicted) whereas a neutral stimulus does not deserve to be recognized. This learning gives emotional "colors" to our world. Some situations are now perceived as associated to a cost (a pain) or a benefit (a pleasure). See [8] [9] for an extensive study of fear conditioning. 3. Then, depending on present internal or external active goals (eg: I am thirsty; I am walking through the room), some stimuli with emotional colors are transformed into motivations (eg: I have to grasp a bottle of beer and it will give me some pleasure; I have to avoid the table or it will give me pain). Using the knoweldge of the consequences of the actions on the world, the correct actions to be performed can be selected. This latter process corresponds to operant conditioning and is clearly related to reinforcement learning. We discuss later how this could lead to less computationally expensive models. The full architecture of our proposal for a control system is shown in figure 1. First, the visual input may trigger direct reflexes after a gross analysis, which represents the innate knowledge of what is absolutely necessary to stay alive. This kind of general architecture is not without common grounds with the Per'Ac architecture introduced in [START_REF] Gaussier | Complex neural architectures for emerging cognitive abilities in an autonomous system[END_REF]. Second, an adaptive visual analysis is carried along two pathways described below.

Visual processing relies on the classical What/Where functional analysis [START_REF] Ungerleider | Object vision and spatial vision: Two cortical pathways[END_REF] [10] : the What pathway represents the analysis of what is currently shown in central vision and the Where pathway associates vision to actions by integrating various informations : the emotional evaluation of what has been seen in central vision, and the position of the object in the visual field (by using the angular position of the eyes). By supposing that the What pathway only deals with central vision, we are simplifying the analysis process and implicitly use ocular saccades to explain how the visual field is analyzed and how to contend with multiple interesting objects. This idea is similar to the embodiment level presented in [START_REF] Dana | Deictic codes for the embodiment of cognition[END_REF] : in the Where pathway, the remaining informations from the What pathway will be the emotional evaluation of the object along with the angular position of the eyes that made the object be in central vision.

Considering the simple scenario where only the color of boxes scattered in the environement matters (blue could mean healthy food while red could be poisoned food), the innate knoweldge would be that when a red box is right in front of the agent, it should avoid eating it and when a box is seen somewhere in the visual field, the agent should get closer. With those two simple reflexes, it is expected to see basic survival behavior : exploration and avoidance of instant death. Over time, the agent will learn that the blue leads to pleasure (satiety) and the red will lead to pain (stomachaches), which is passive conditioning. It will also learn that when there is something that leads to pleasure on the right, it should turn right, while if it is something painful, it should go to the left to avoid it, which is operant conditioning.

As a proof of concept, we now present how each of the modules of this system has been implemented.

Visual analysis

It is not our purpose to propose here novel methods of analysis of the visual input. The Where pathway in our model is not purely about location, it is also a How pathway considering how it leads directly to action. Our What analysis extracts shape (orientation in our experiments) and color, as in V4. Our Where pathway extracts the location of the object currently seen in central vision in an egocentric coordinate system, using among other informations the angular position of the eyes.

Classification modules

The What and Where classifications use the exact same implementation. The main purpose of this system is, as shown in figure 2, to classify an arbitrarily large input vector into as few cases as possible. It is an incremental neural network similar to Fritzke's Growing Neural Gas as described in [START_REF] Fritzke | A growing neural gas network learns topologies[END_REF] without the topology as it is irrelevant in this context. We decided to use incremental learning as it is clear that it is not possible to consider all possible configuration of stimuli. For instance, in the What pathway, the characteristics extracted from the visual input are mainly shapes and colors : the number of neurons necessary to encode each possible configuration will increase heavily as more colors or shapes are considered, even though not all configurations are relevant. This is why we have added two inputs to the classification modules : novelty and importance. The novelty of a configuration in the case of the What pathway is computed from the difference between the current emotional evaluation and the real pain and pleasure signals. In other words, if the body is feeling pleasure from eating while the emotional evaluation system is not predicting it, the novelty signal will be high. Likewise, if the emotional evaluation is predicting pain while the body is feeling pleasure, the novelty will also be high. The novelty signal is irrelevant for the Where pathway, it has been set to a high value to let the body react more quickly to objects it has already identified. In both pathways, the importance signal is the sum of the pain and pleasure signals coming from the body : if it is high, it means that something is going on and learning should take place. Those two signals are used to control the creation rate of new neurons to represent configurations of input stimuli and to control how close the prototype of the neurons will be to the input stimuli. Unlike the vigilance signal in Grossberg's ART [START_REF] Carpenter | Adaptive resonance theory[END_REF], which shares many common points with our classification modules, our novelty and importance signals are tightly bound to the rest of the system and are not external parameters chosen in advance.

Conditioning modules

The pain and pleasure conditioning modules are based on Moren's amygdala model, as presented in [START_REF] Jan | Emotion and Learning : a computational model of the Amygdala[END_REF] , extended to allow real-time learning by adding a Short Term Memory trace to the inputs. As explained previously, the input of the conditioning modules are not directly the visual stimuli but a small subset of configurations chosen by the classification module. Concerning the non-stationary nature of the classification, we also had to modify Moren's model to unlearn association which have lost their meaning due to the modification of the classification. Specifically, when the classes in the what classification module change, the association of the inputs in the conditioning modules with the emotional response decrease.

Winner Takes All

Each of the output neuron of the classification in the where pathway is linked with a weight between 0 and 1 to each neuron in the Winner take all module. The neuron of the WTA which has the highest output becomes the only one firing to control the associated action. To learn those weights, we use a mechanism very similar to Q-Learning using TD(0). Once a neuron of the ouput layer of the classification has been activated, it has a very slow decay, hence keeping a trace of the fact that it has been chosen for a long time. In the TD(0) algorithm, this neuron would encode the eligibility of the action. When pleasure or pain is felt, weights increase or decrease accordingly, hence favoring or leading to avoidance of the last chosen actions, proportional to the value of the eligibility trace.

On a side note, the reflex pathway is connected to the WTA with small constant weights so that a learned behavior can quickly but not instantaneously take over reflexes.

Discussion

The architecture has been simulated and is now being tested on a wheeled robot, with the scenario given in the architecture overview. We expect to observe the mentioned bootstrapping effect and autonomous learning once we have adapted all the time constants and various parameters from the simulation to the actual robot.

This architecture offers a base on which higher levels of abstraction can build upon. It is an attempt at creating a comprehensive learning and control system based on biological data and vastly supported hypothesis. In a similar way than how the What/Where pathways are built upon the Reflexes pathway, we strongly believe that it will be possible to implement higher level functionalities which will take over those reflexes when necessary, such as the executive functions of the frontal lobe [3] [5].

It is interesting to note how abstract concepts such as frustration or relief (difference between the expected pleasure/pain and the actual pleasure/pain, called importance in the classification module), fear and anticipation quickly stem from the embodiment of the learning process. Our point is not to decide whether those are actual feelings or not, but how the use of the semantics of emotions can be helpful in designing learning system in a bottom up approach.

As it is a work in progress, many improvements are expected. Currently, our main focus is on the Where pathway : as mentioned earlier, it is our belief that one of the ways to deal with the computational complexity, that usually comes with control systems learning by classical reinforcement learning, is to firstly only consider pertinent elements (emotionaly colored) in the environment, and secondly only consider the actions related to the current movitations. Such ways of computing are possible, and even natural, with an embodied cognition and richer internal states.

Fig. 1 .

 1 Fig. 1. Overview of the architecture

Fig. 2 .

 2 Fig. 2. Schematic representation of the classification modules

Corresponding author: frederic.alexandre@loria.fr


Ichiro Imanishi

Aurélie Nicolas

Ana-Carolina Barbosa Caetano

Thiago Luiz

Paula Castro

Natayme Rocha Tartaglia

Ricardo Mariutti

Eric Guédon

Sergine Even

Nadia Berkova

Raghuvir K Arni

exfoliative toxin e, a new Staphylococcus aureus virulence factor with host-specific activity

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

In the Supplementary Information file originally published with this Article, supplementary figure S6, showing Best HADDOCK solutions for ETE docking on Dsg1 of four animal species, and supplementary table S4, showing HADDOCK scores, were omitted. These errors have been corrected in the Supplementary Information that now accompanies the Article.

Due to the omitted figure and table, there were additional errors in the main body of the Article.

In the Results section, under the subheading 'Two types of docking orientation predicted' , "Docking simulations were also carried out with human and canine desmogleins (Supplementary Fig. S1) and displayed the same two ETE orientations. " now reads, "Docking simulations were also carried out with human and canine desmogleins (Supplementary Fig. S6) and displayed the same two ETE orientations. "

In the Experimental Procedures section, under the subheading 'Recombinant ete gene product' , "Purity of the ETE protein was determined by SDS-PAGE gels (Fig. S2). " now reads, "Purity of the ETE protein was determined by SDS-PAGE gels (Fig. S1). "

Under the subheading 'In vitro digestion of Dsg1 by ETs' , "Non-cropped, non-modified immunoblotting images for in vitro digestion of recombinant Dsg1 are presented in supplemental data (Fig. S3). " now reads, "Non-cropped, non-modified immunoblotting images for in vitro digestion of recombinant Dsg1 are presented in supplemental data (Fig. S2). "

Published: xx xx xxxx
open "Residues involved in the active sites of both proteins were required. The active ETE residues were those defined in 9 (His96, Asp145 and Ser219) and those for Dsg1 were those defined in 52 (Glu381, Gly382 in protein sequence matching Glu332, Gly333 in structural models). " now reads, "Residues involved in the active sites of both proteins were required. The active ETE residues were those defined in Mariutti et al. 9 (His96, Asp145 and Ser219) and those for Dsg1 were those defined in Hanakawa et al. 52 (Glu381, Gly382 in protein sequence matching Glu332, Gly333 in structural models). "

  Staphylococcus aureus virulence factor with host-specific activity ichiro imanishi, Aurélie nicolas, Ana-carolina Barbosa caetano, thiago Luiz de paula castro, natayme Rocha tartaglia, Ricardo Mariutti, eric Guédon , Sergine even, nadia Berkova , Raghuvir K. Arni, nubia Seyffert, Vasco Azevedo, Koji nishifuji & Yves Le Loir Correction to: Scientific Reports https://doi.org/10.1038/s41598-019-52777-3, published online 08 November 2019

Scientific RepoRtS | (2020) 10:2460 | https://doi.org/10.1038/s41598-020-59360-1 www.nature.com/scientificreports

Scientific RepoRtS |(2020) 10:2460 | https://doi.org/10.1038/s41598-020-59360-1