Sebastián Nash

Cristián Escauriaza

Pascal Finaud-Guyot

Wolfram Jahn

Antoine Rousseau

Single Porosity Model: Exploring the Spatial Resolution Limits in Complex Urban Patterns

Keywords: urban floods, porosity shallow water equations, flood hazard, numerical modeling, Godunov-type schemes

When modeling large-scale urban floods, the use of porosity non-linear shallow water equations (PNSWE) emerges as an interesting sub-grid approach for reducing computation time while preserving the structure of the solution. In such models, fine-scale topographic information is represented at a coarser scale through porosity parameters, enabling a speed-up in computations at the expense of losing accuracy while computing hydrodynamic variables. Porosity accounts for both the change in storage and exchange section due to the presence of obstacles (e.g., houses or buildings) in the floodplain and introduces a new source term associated with its gradient into the momentum equations. Here, we use the Single Porosity model (SP) in Cartesian coordinates to simulate flows in both an idealized and a real-world urban area, while gradually increasing the spatial resolution. During such partial coarsening, in which we move from fine-to macro-scale, the porosity distribution changes within the urban zone from a highly heterogeneous field

Introduction

Floods stand as the world's most recurrent disaster, inflicting significant damage, casualties, and economic losses every year [START_REF] Liu | Association between floods and typhoid fever in yongzhou, china: Effects and vulnerable groups[END_REF]. Recent investigations have shown that between 1980 and 2020, more than 4800 flood disasters have occurred, resulting in an annual average of more than 5900 deaths and $21 billion USD in losses [START_REF] Duan | Assessment and spatiotemporal analysis of global flood vulnerability in 2005-2020[END_REF]. These events have also been responsible for almost half of the human fatalities from natural hazards [START_REF] Lee | Assessing the vulnerability of different age groups regarding flood fatalities: case study in the Philippines[END_REF][START_REF] Fitzgerald | Flood fatalities in contemporary australia (1997-2008)[END_REF].

Analyses of future flood risks and associated damages show that they will be strongly influenced by two critical factors: urbanization and climate change [START_REF] Barría | From multi-risk evaluation to resilience planning: The case of central chilean coastal cities[END_REF][START_REF] Pariartha | Compounding effects of urban-ization, climate change and sea-level rise on monetary projections of flood damage[END_REF][START_REF] Wang | Assessing urban flooding risk in response to climate change and urbanization based on shared socio-economic pathways[END_REF]. As urbanization continues to accelerate, an increasing settlement of people in high-risk prone areas has led to a higher exposure of population and in-frastructure to potential flooding events [START_REF] Yari | Underlying factors affecting death due to flood in iran: A qualitative content analysis[END_REF][START_REF] Alderman | Floods and human health: A systematic review[END_REF][START_REF] Jonkman | Loss of life due to floods[END_REF]. This has been exacerbated by improper land-use planning [START_REF] Jha | Cities and flooding: a guide to integrated urban flood risk management for the 21st century[END_REF], rapid transitions in land-use that limit or prevent water infiltration in cities [START_REF] Jodar-Abellan | Impact of land use changes on flash flood prediction using a sub-daily swat model in five mediterranean ungauged watersheds (se spain)[END_REF][START_REF] Lee | Development of a flood-damage-based flood forecasting technique[END_REF], and critical infrastructure near river channels exposed to scour in subcritical or supercritical conditions [START_REF] Link | Scour at bridge foundations in supercritical flows: An analysis of knowledge gaps[END_REF][START_REF] Nieto | Fragility curves for road embankments exposed to adjacent debris flow[END_REF].

In addition, climate change-driven alterations in rainfall patterns and sea level rise, will affect the magnitude, frequency, and intensity of urban flood events, which are expected to increase continuously (IPCC, 2021;[START_REF] Jamali | A rapid urban flood inundation and damage assessment model[END_REF][START_REF] Chen | Urbanization and climate change impacts on future flood risk in the pearl river delta under shared socioeconomic pathways[END_REF]. Climate modifications are also expected to influence interannual variations associated to climatic oscillations, such as the El Niño Southern Oscillation (ENSO), which may induce strong spatiotemporal fluctuations in flood occurrences [START_REF] Ward | Annual flood sensitivities to el niño-southern oscillation at the global scale[END_REF].

Flood risk assessment (FRA) is a critical feature for flood risk management and it plays an essential role in informing decision-making and the formulation of effective mitigation measures [START_REF] Zou | Comprehensive flood risk assessment based on set pair analysis-variable fuzzy sets model and fuzzy ahp[END_REF]. Effective FRA approaches must be based on reliable estimates of flood hazards, exposure, and vulnerability. An important element of the first two components corresponds to the computation of the flow properties within urban areas, which define the magnitude of the event on-site [START_REF] Bruwier | Shallow-water models with anisotropic porosity and merging for flood modelling on cartesian grids[END_REF]. Therefore, the requirement to assess flood hazards needs precise modeling tools to evaluate the propagation of floods in urban systems [START_REF] Ferrari | Floodwater pathways in urban areas: A method to compute porosity fields for anisotropic subgrid models in differential form[END_REF], making numerical models a vital component of FRA [START_REF] Rosenzweig | The value of urban flood modeling[END_REF][START_REF] Sanders | Integral formulation of shallow-water equations with anisotropic porosity for urban flood modeling[END_REF].

The classical Nonlinear Shallow Water Equations (NSWE) or Saint-Venant equations have been generally considered state-of-the-art for the simulation of inundations and flood paths, replacing one-dimensional former approaches [START_REF] Costabile | Enhancing river model set-up for 2-d dynamic flood modelling[END_REF]. A main concern with such models is the extensive computational time required to numerically solve the system in difficult topographies with the necessary resolution in extensive urban areas, which can be in the order of centimeters [START_REF] Bates | Integrating remote sensing data with flood inundation models: how far have we got?[END_REF][START_REF] Özgen | Shallow water equations with depth-dependent anisotropic porosity for subgrid-scale topography[END_REF][START_REF] Kim | Urban flood modeling with porous shallow-water equations: A case study of model errors in the presence of anisotropic porosity[END_REF]Contreras and Escauriaza, 2020).

Performing a simple coarsening of structured meshes generally leads to large errors in flood predictions because it degrades the representation of buildings and other infrastructure [START_REF] Yu | Urban fluvial flood modelling using a two-dimensional diffusion-wave treatment, part 1: mesh resolution effects[END_REF][START_REF] Brown | Modeling storm surge flooding of an urban area with particular reference to modeling uncertainties: A case study of canvey island, united kingdom[END_REF].

The fundamental issue is that computational meshes can only resolve features considerably larger than the grid size [START_REF] Sanders | Integral formulation of shallow-water equations with anisotropic porosity for urban flood modeling[END_REF]. Thus, Porosity Nonlinear Shallow Water Equations (PNSWE) emerge as an interesting and consistent sub-grid approach for modeling urban floods, where the porosity accounts for both the change of the storage and exchange section due to the presence of buildings in the inundated areas [START_REF] Dewals | Porosity models for large-scale urban flood modelling: A review[END_REF].

In models based on the PNSWE, fine-scale topographic information is represented at a coarser scale using porosity parameters, allowing a significant decrease in computational time at the expense of losing accuracy in the calculation of hydrodynamic variables. The early approaches to PNSWE emerged in a differential formulation of the equations, like micro-topography sub-grid modeling [START_REF] Defina | New set of equations for very shallow water and partially dry areas suitable to 2d numerical models[END_REF]. These early models inspired later developments focused on urban areas, such as the Single Porosity (SP) model [START_REF] Guinot | Flux and source term discretization in two-dimensional shallow water models with porosity on unstructured grids[END_REF]. Building upon this, further advancements addressed urban terrain anisotropies, as seen in models such as the Multiple Porosity (MP) model [START_REF] Guinot | Multiple porosity shallow water models for macroscopic modelling of urban floods[END_REF] and the addition of a building drag formulation [START_REF] Velickovic | Steady-flow experiments in urban areas and anisotropic porosity model[END_REF]. Simultaneously, PNSWE models were developed based on integral formulations, primarily focused on representing porosity using local geometric parameters. The pioneering work in this direction is the Integral Porosity (IP) model developed by [START_REF] Sanders | Integral formulation of shallow-water equations with anisotropic porosity for urban flood modeling[END_REF]. This model provided the foundation for the development of subsequent models, including the Dual Integral Porosity (DIP) model [START_REF] Guinot | Dual integral porosity shallow water model for urban flood modelling[END_REF], which aimed to overcome mesh-dependency issues associated with the IP model. Moreover, advancements such as incorporating depthdependent porosity, as demonstrated in the Depth-Dependent Dual Integral Porosity (DD-DIP) model of [START_REF] Guinot | Flux closures and source term models for shallow water models with depth-dependent integral porosity[END_REF], further enriched the integral approach.

The SP model, based on the differential form of the equations, assigns a unique porosity value per cell for large-scale resolutions to maintain the parsimony of the formulation. This porosity, referred to as storage porosity (ϕ), corresponds to the fraction of the plan view area available to the flow, while (1 -ϕ) represents the fraction of the surface occupied by buildings. The SP model is particularly interesting since it retains mathematical properties of the classical NSWE model such as hyperbolicity, rotational invariance, and celerity of small disturbances [START_REF] Varra | On integral and differential porosity models for urban flooding simulation[END_REF]. The first SP models assumed a single constant value of the porosity parameter for the cells within the urban area. The porosity is then computed at a large-scale level, considering all the extent of the city, and then imposed on the cells within the urban domain, while the outer cells have an assigned value of ϕ = 1. As a consequence, the source term involving the porosity parameter only accounts for a local head loss at the entrance and exit of the urban area. The head losses occurring inside urban areas are commonly accounted for by using drag force-type of source terms [START_REF] Soares-Frazao | Investigation of distributed-porosity fields for urban flood modelling using single-porosity models[END_REF].

As pointed out by [START_REF] Sanders | Integral formulation of shallow-water equations with anisotropic porosity for urban flood modeling[END_REF], the formulation of differential equations describing flow in porous media is accomplished by postulating the existence of a Representative Elementary Volume (REV), which is required for continuum purposes. The REV, now understood as a representative area rather than volume (used for consistency with previous published works), refers to an area large enough to capture average properties that can be represented as porosity in this case and thus characterizes the entire urban domain. Therefore, this uniform porosity value should be computed at the REV level. Initial work on the definition of REV demonstrated that for typical European capital cities, a large grid size of at least 2.3 km is needed [START_REF] Guinot | Multiple porosity shallow water models for macroscopic modelling of urban floods[END_REF]. In many situations, however, this size may exceed the spatial scale of the urban area itself, making the REV meaningless for small towns or urban areas that are spatially limited. Practical applications of the PNSWE, however, have demonstrated the viability of finite-volume schemes for the solution of the SP model with a local storage porosity defined at the level of computational cells [START_REF] Velickovic | Shallowwater model with porosity: sensitivity analysis to head losses and porosity distribution[END_REF][START_REF] Özgen | Wave propagation speeds and source term influences in single and integral porosity shallow water equations[END_REF][START_REF] Soares-Frazao | Investigation of distributed-porosity fields for urban flood modelling using single-porosity models[END_REF]. In the same line, and as discussed by [START_REF] Guinot | Multiple porosity shallow water models for macroscopic modelling of urban floods[END_REF], the existence of the REV in the urban area has no consequences with respect to the applicability of the equations.

Numerical studies using porosity at the computational cell-level, which implies defining a non-uniform porosity distribution within the city, have im-proved the results compared to the constant porosity value inside the urban domain [START_REF] Soares-Frazao | Investigation of distributed-porosity fields for urban flood modelling using single-porosity models[END_REF]. In some approaches, the original porosity definition, applied to the cell-level, has been manipulated further in order to represent the city pattern and improve the flow depth estimation, like the building-peak distribution proposed by [START_REF] Linkens | Application du concept de porosité pour la modélisation d'écoulements en milieu urbain[END_REF].

In this investigation we use the SP model within a Cartesian frame to solve the propagation of urban floods, employing the original storage porosity (ϕ) definition at the computational cell-level, and evaluating the consequences of varying the size of the discretization. Motivated by the development of earlywarning systems and faster-response models to analyze multiple potential scenarios before a extreme hydrometeorological event, the primary objective is to evaluate predictions of water depths, flood limits, and hazard index for various spatial resolutions, providing new insights on the advantages and limitations of this formulation.

It is also of special interest to understand the trade-offs between the reduction of the computational time and the prediction error induced by the coarsening process. From this systematic procedure, we study the spatial resolution limits for the SP model to deliver reliable results that can be trusted by decision makers and managers. This formulation can be implemented in early-warning systems that typically require a rapid response and the evaluation of multiple potential scenarios, based on a physical framework and without relying only on machine learning approaches.

The model is applied in two cases: (i) An idealized non-periodic city; and

(ii) A real-world case study located in La Reina, a district within the city of Santiago de Chile. A fine-scale well-balanced, widely validated NSWE model, introduced by Guerra et al. ( 2014) is used for computing a reference solution (see also [START_REF] Loose | Finite volume of variable shallowwater flow equations for well-mixed estuary: application to the rio Maipo estuary in central Chile[END_REF]. These high-resolution calculations provide detailed results (Contreras and Escauriaza, 2020), but they are highly time-consuming at the urban resolution obtained from LiDAR or satellite images. In our analysis, we systematically analyze the effects of downgrading the spatial resolution, incorporating the sub-grid effects through the PNSWE model.

This paper is organized as follows: In Section 2, we present the approach for the numerical modeling of urban flood, detailing the governing equations, the collection of numerical methods used in the calculation, and the mesh coarsening process. In Section 3, we describe the two different simulation cases and the mesh resolutions employed in each case. Section 4 provides numerical results and discussion of computational time, incorporating the analysis of the physical variables that define the event, such as the maximum depth, flood limits, and hazard estimation. Finally, in Section 5, we present the conclusions and closing remarks.

Modeling Urban Floods with Porosity

Governing Equations

The shallow water equations correspond to a simplification of the threedimensional Navier-Stokes equations, where the main assumptions are: (i) incompressible and homogeneous fluid; (ii) hydrostatic pressure distribution;

(iii) negligible vertical velocities; and (iv) vertically averaged horizontal velocities. When adding storage porosity to the classical model and using the differential form, we obtain the PNSWE. In matrix form, they are expressed as follows:

∂ ∂t (ϕQ) + ∂ ∂x (ϕF ) + ∂ ∂y (ϕG) = S B (Q) + S F (Q) + S P (Q), (1)
with

Q =      h hu hv      ; F (Q) =      hu hu 2 + gh 2 2 huv      ; G(Q) =      hv huv hv 2 + gh 2 2      ; S B (Q) =      0 -gh ∂z b ∂x -gh ∂z b ∂y      ; S F (Q) =      0 -S f x -S f y      ; S P (Q) =      0 gh 2 2 ∂ϕ ∂x gh 2 2 ∂ϕ ∂y      , ( 2
)
where ϕ is the storage porosity, h corresponds to the water depth, u and v are the depth-averaged flow velocities in the x and ŷ Cartesian directions, respectively, t indicates time, z b is the bed elevation, g is the gravity acceleration (9.81m/s 2 ) and S f x and S f y correspond to the friction source term in the x and ŷ direction, respectively. Note that a distinction is made between the source terms, where S B (Q) accounts for the bed slope, S F (Q) is the bed and wall friction, and S P (Q) is the porosity gradient, which appears when adding the porosity term to the classical model.

Numerical Model

The numerical solution of the system of equations ( 1) is obtained by implementing a second-order Godunov-type scheme, where time stepping is performed by means of a four-stage Runge-Kutta scheme [START_REF] Guerra | Modeling rapid flood propagation over natural terrains using a wellbalanced scheme[END_REF]. For the resolution of the Riemann problem, we use the original SP formulation, which follows a lateralized discretization for the pressure source terms in the momentum flux. Therefore, the static equilibrium conditions are satisfied even when the porosity is different on both sides of the interface [START_REF] Soares-Frazao | Investigation of distributed-porosity fields for urban flood modelling using single-porosity models[END_REF]. The source term that accounts for friction effects is discretized using a time-splitting procedure. Thus, the state variable updating involves two stages: one for considering the friction source term and the other for addressing the hyperbolic effects.

On the first step, the semi-discrete formulation of the system results in:

d dt (ϕQ) = S F (Q). ( 3
)
Using Manning equation to parameterize flow resistance, and linearizing the system with respect to the known time level n, leads to a system of ordinary differential equations with an analytical solution [START_REF] Soares-Frazao | Investigation of distributed-porosity fields for urban flood modelling using single-porosity models[END_REF].

With these updated velocities, the friction-modified variables are computed, which serve as input for the next time step.

On the second step, the system to solve considers only the bottom slope and porosity source term, leading to:

∂ ∂t (ϕQ) + ∂ ∂x (ϕF ) + ∂ ∂y (ϕG) = S B (Q) + S P (Q). (4)
The finite-volume semi-discrete system can be written as follows,

d dt (Q i,j ) = 1 ϕ i,j ∆x i,j (ϕF i-1/2R,j -ϕF i+1/2L,j )+ 1 ϕ i,j ∆y i,j (ϕG i,j-1/2R -ϕG i,j+1/2L ), (5)
where ϕ i,j and Q i,j are cell-centered values, ∆x i,j and ∆y i,j are the computational cell widths in the x and ŷ directions, respectively, and ϕF and ϕG are the numerical fluxes obtained by solving the RP in the x and ŷ directions, respectively. The sub-index i-1/2 denotes the cell's left-side interface, i+1/2 the right-side interface, j -1/2 the lower-side interface, and j +1/2 the upper-side interface. The subscript L indicates the left-side of the interface (lower-side in ŷ), while R indicates the right-side (upper-side in ŷ).

In the initial part of this step, and with the aim of achieving a secondorder accurate discretization, a MUSCL reconstruction is performed for the primitive variables h, u, v, and η, the last one being the free-surface elevation (η = z b + h), using the minmod slope limiter. Subsequently, the RP is solved with an HLL-modified solver, which accounts for a specific treatment for the continuity and momentum equations. The source terms S B (Q) and S P (Q) discretizations are included in the Riemann solver; therefore, their effect is accounted for through the fluxes ϕF and ϕG.

The two-dimensional CFL condition is used to calculate the time step while ensuring numerical stability, the same stability criterion used for the classical models, as the celerity of the waves is not affected by the inclusion of a single porosity value per cell. The time step is computed as:

∆t n = min(∆x i,j , ∆y i,j ) max(max(u n i,j + gh n i,j , v n i,j + gh n i,j ))

• CFL. (6)

Meshing: From fine-to macro-scale

When modeling urban floods, we distinguish among three different modeling spatial scales [START_REF] Bruwier | Shallow-water models with anisotropic porosity and merging for flood modelling on cartesian grids[END_REF]):

• Fine-scale: all obstacles in urban areas are resolved explicitly and with a fairly good accuracy (about ten cells over the typical width of flow paths).

• Meso-scale: the obstacles are explicitly represented in the computational domain, but poorly resolved (just a couple of cells over the typical width of flow paths).

• Macro-scale: the obstacles have a typical size comparable or smaller than the mesh size (flow paths are generally smaller than the size of a single computational cell).

For the reference solution computed with the classical NSWE, we use finescale modeling. In an initial step, we treat the bare-earth topography and the buildings separately. For the bare-earth topography, we use a Digital

Terrain Model (DTM), which is commonly obtained by means of LiDAR measurements. To account for the effects of the buildings on the flow, we use the Building-Block (BB) representation, where the bed elevation is artificially increased in the area where a building is situated. By setting an artificially high elevation, large enough to ensure that no water overflows the buildings, the flood would only flow around them. Such information about the location of the buildings is provided by the Binary Buildings Map (BBM). The BBM corresponds to the field built from the discrete binary variable B(i, j), where B(i, j) = 1 indicates a void or street, while B i,j = 0 denotes a building. The geometry of the BB domain is therefore obtained by combining the BBM and the DTM.

For the SP model meshing, the structured fine-scale mesh is coarsened by merging neighboring cells (in an equal number of rows and columns). In such a way, the porosity at the computational cell-level for a coarse cell (i, j)

that contains W fine cells merged is computed from the BBM as

ϕ i,j = w∈W [(B) w • (A) w ] w∈W (A) w , (7)
where A w corresponds to the area of the wth merged cell. Note that the original definition of storage porosity ϕ is preserved at the cell-level. Similarly to porosity, the resulting coarse terrain elevation is computed as:

(z b ) i,j = w∈W [(z b ) w • (A) w ] w∈W (A) w . ( 8
)
Since we are working within a Cartesian reference frame and from a practical point of view, topography coarsening consists of DTM space averaging, while porosity computation consists of BBM space averaging.

Simulations Cases

In this section, we briefly summarize the two cases analyzed in this investigation. We first study the performance of the SP model in an idealized small domain, which is designed with a constant slope and a non-symmetric heterogeneous distribution of buildings. In this first case the streets are To evaluate the model performance in this case, we carry out simulations considering a dam-break flow introduced in the direction of the arrow shown in Figure 1, which is the most severe condition that can be tested in these urban-like domains. For the remaining downstream boundaries, a free-exit condition is defined [START_REF] Guerra | Modeling rapid flood propagation over natural terrains using a wellbalanced scheme[END_REF]. As an initial condition, the dam has 1m of water depth, while a dry-bed is used for the rest of the domain. The fine-scale mesh has a resolution of 1m in order to have at least 9 cells over the typical width of the flow paths, fitting the fine-scale definition.

For the partial transition to the macro-scale, we use the following resolutions (in meters): 2, 3, 4, 6, and 12. In the case of res=12m where the porosity map is shown in Figure 1b), the cell size is larger than SW, thus fitting the macro-scale definition.

Case 2: La Reina, Santiago de Chile

La Reina is a district situated in the Metropolitan area of Santiago, Chile, on the eastern border of the city. On this edge, and as shown in Figure 2, it limits with the Andean foothills, specifically with the Quebrada de Ramón basin, in which the main channel flows directly into the urban district (see Contreras and Escauriaza, 2020, for details). During May 1993, one of the most significant flood event recorded in this area occurred, triggered by a warm storm associated with the ENSO phenomenon [START_REF] Garreaud | Análisis meteorológico de los aluviones de antofagasta y santiago de chile en el periodo 1991-1993[END_REF], that elevated the 0 • C isotherm up to 4000m masl (Contreras and Escauriaza, 2020).

These conditions generated high-altitude heavy rainfall in the Andes mountains, and as a consequence, a substantial debris flow rapidly inundated the district in a matter of minutes. The event caused the unfortunate loss of 26 lives, 8 missing people, and major infrastructural damage, with a total financial cost of $5 USD million. The reported affected area is shown in Figure 2 [START_REF] Naranjo | Flujos de detritos y barro que afectaron el sector oriente de Santiago, el 3 de mayo de[END_REF], denoted by a yellow marking.

We perform numerical simulations in this urban section, with a computational domain that covers the green area of the urban zone, as shown also in Figure 2. For the topography, a 1m resolution DTM of the area is used, obtained through LiDAR sensing. Using data from OpenStreetMap (OSM), the BBM is built by subtracting the streets and park areas from the resi- dential area. With this procedure, and from a practical perspective, we are assuming that everything that is not a street or park is a building. Thus, the characteristic SW corresponds to 14m for common streets, while for main streets, SW corresponds to approximately 30m.

The DTM indicates a mean slope of around -1°in the E-W direction and -0.65 • in the N-S direction. A dam-break flow is introduced in the area where the channelized stream overflows, according to previous simulations from Contreras et al. (2020), following the direction indicated in Figure 3.

For the remaining downstream boundaries, a free-exit condition is defined.

As initial condition, the dam also has 1m of water depth, while a dry-bed is used for the rest of the domain. Since the refined mesh resolution is 1m, we have at least 14 cells over the width of the flowpaths, adjusting to the fine-scale definition. For the partial transition to the macro-scale, we use the following resolutions (in meters): 4, 8, 12, 32, and 64. The porosity field for res=32m is shown in Figure 3b).

Results and Discussion

Computational Time

Here we compare the CPU time of the reference model with the SP model at different spatial resolutions to quantify the computational advantage of SP.

For case 1, shown in Table 1, a significant reduction is achieved when using the SP model, which is even more considerable for higher spatial resolutions.

The percentages of CPU time, with respect to the reference model, vary between 12% for a resolution of 2m and 0.1% for a resolution of 12m.

For case 2, where the CPU times are shown in Table 2, there is also a notable reduction in computation time. When using a resolution of 4m, the CPU time is around 16% with respect to the reference model. This reduction becomes more noticeable when using higher resolutions, reaching a CPU time of 0.04% for the 64m resolution.

We can understand the reduction in time when using the SP model as a combined effect of two factors. First, by using higher spatial resolutions for the same domain, we have fewer cells and, therefore, fewer interfaces of the discretization where the Riemann problem has to be solved. This effect always yields a time reduction when increasing spatial resolution. In addition, the time-step is less constrained by the cell dimensions according to the stability condition stated in equation ( 6) as the area per cell is larger.

However, it is important to keep in mind that when using coarse grids with the SP model, there is a tendency to overestimate velocities within the urban domain. This generates the opposite effect, decreasing the time-step. Thus, it is not possible to know a priori the effect of mesh resolution over the time-step.

For case 1, we observe a clear tendency to increase the mean ∆t with spatial resolution. In this way, both effects act in the same direction to reduce the total computational time.

For case 2, we note that, even when all values of the mean ∆t of SP are greater than the one from the reference model, there is not a clear tendency with respect to the mesh resolution. For example, when comparing the 32m with the 64m resolution, we notice that the former has a mean ∆t around 4 times larger, but still presents a reduced total time with respect to the latter. This suggests that the effect of computing less Riemmann problems at interfaces is more relevant.

Finally, it is important to consider that the time reduction achieved through increased spatial resolution must be evaluated along with the error induced by the SP model. Although further time reduction is achievable with coarser meshes, it is crucial to be aware that the hydrodynamic variables lose physical meaning, as we discuss in the following subsection.

Maximum water depth

Here we compute the maximum water depth field through the total simulation time for each spatial resolution in the streets. With the objective of having an overview of the mean induced error, we compute the Mean Absolute Error (MAE). Also, we consider the Mean Bias Error (MBE) in order to assess the systematic bias of the depth estimation in terms of underestimation or overestimation. In Figure 4, we observe a general trend of increasing MAE with the mesh resolution that can be anticipated due to the topographic detail degradation.

The MBE indicates a systematic underestimation of the maximum water depth in the streets. This suggests that the SP model has a more diffusive character when the mesh size increases. Such results can be explained since in the SP model obstacles lose their representation as the mesh size increases, the water spreads over a larger area and converges to a lower water depth. As a consequence, when considering only the flow in the streets, a lower water depth is found.

In spite of these apparent shortcomings, the MAE remains below 4cm for the 12m resolution in case 1, while for case 2, it remains under 5cm when using a 64m mesh. The systematic underestimation in the streets is always less than 1cm on average for case 1. For case 2, the underestimation increases with mesh resolution, presenting an average underestimation of up to 3cm for a 64m resolution mesh. In both cases, the errors presented are on the order of a few centimeters for the maximum water depth, which are remarkable results in the context of a drastic change in spatial resolution.

It is important to consider not only the numerical averages but also the spatial distribution of the flood within the city. For instance, we could think that using a 64m grid and having an average underestimation of 3cm is an encouraging result for La Reina in case 2, even more so when considering the reduced computational time of the SP model option. However, as we will analyze in the following section, the model may lose its physical coherence when a very coarse resolution is employed.

Flood Extension

With the objective of analyzing the spatial distribution of the flood, we compute the inundated region limits for a threshold of 0.1m. For the case 1, the flood map shown in Figure 5, demonstrates that the inundation pattern is captured, especially up to a resolution of 6m. Note that the flood pattern from the reference is represented for all SP resolutions below street width, which is 9m in this case.

Here, the effect of the resolution on the development of a shallower and more extended flow is well illustrated (e.g., while comparing 5a) with 5f)),

For intermediate cases, and specially for the 4m and 8m resolutions, where the porosity gradient becomes more relevant, it can be seen that such directionality is partially recovered. These results suggest that, similarly to case 1, for resolution values below SW (14m in this case), the flood patterns can be successfully captured in general terms.

A more detailed view of the mentioned north section is shown in Figure 7, where the color green indicates a prediction agreement with the refined case, red corresponds to an underprediction from SP, and blue an overprediction.

In this analysis is clear that for the cases with resolution below SW, the SP model systematically fails to capture the inundation area. This seems to be in agreement with the fact that the mesh manages to resolve features that are larger compared to the cell size.

Hazard Index

For conducting a correct FRA, in addtion to a thorough vulnerability analysis, it is crucial to identify and correctly estimate the current and possible areas exposed to flood hazards [START_REF] Erena | Flood hazard mapping using flo-2d and local management strategies of dire dawa city, ethiopia[END_REF]. Thus, flood hazard maps can be used to effectively reduce associated damage [START_REF] Ogato | Geographic information system (gis)-based multicriteria analysis of flooding hazard and risk in ambo town and its watershed, west shoa zone, oromia regional state, ethiopia[END_REF]. For computing flood hazards, we use criteria based on the maximum water depth and flow velocity, as proposed by [START_REF] Garcia | Debris flows of december 1999 in venezuela[END_REF].

q = h √ u 2 + v 2 as: (i) Low for 0.1m < h < 0.5m ∧ 0.1m 2 /s < q < 0.5m 2 /s;

(ii) Medium for 0.5m < h < 1.5m ∨ 0.5m 2 /s < q < 1.5m 2 /s, and (iii) High for h > 1.5m ∨ q > 1.5m 2 /s. The hazard map for case 1 in Figure 8 shows that the main structure of the Low hazard area is well captured up to a 6m resolution. We can see that the plumes leaving the urban area from the sides are also captured at this resolution. For the case of 12m resolution, there is a significant overestimation of the flooded area in the urban zone surroundings, especially in the lateral boundaries near the downstream exit. For the Medium hazard area, we can see in the reference model three main branches entering the city along the main flow-aligned streets and one deviated branch to the city front limit. Such branches are captured at all scales, where the 12m model fails to predict some minor areas inside the urban domain. For the High hazard areas, in red, we identify three locations from the reference: in the flow inlet, in the section of the urban area directly exposed to the flow, and in the center at a dead-end street. Only the section facing the dam-break flow is consistently captured for all resolutions, while the other zones are captured up to a 4m resolution. Overall, the main structure of the hazard zones and their spatial distribution are in good agreement with the high-resolution reference case up to a 6m resolution.

For case 2, shown in Figure 9, we observe a clear tendency to underestimate the hazard by the SP model. When using a 4m or 8m resolution, the general location of hazard zones is captured, but many streets in the Medium hazard category are estimated as Low hazard with SP. For narrower streets, a High hazard level is found in the refined model because of increased flow velocities, which are captured with SP until a 16m resolution. However, for resolutions from 16m and beyond, as preferential flow paths are not captured, there are many zones with no hazard estimated. Thus, the inability of SP to accurately capture these paths directly affects the hazard assessment, evidencing major discrepancies with the refined model. Here it is very clear, observing the 32m and 64m resolutions, that the flood advances conforming to the bed slope as the resolution is increased, making the hazard map lose its physical meaning. (2005). The systematic coarsening shows significant differences on the classification as the grid size is increased.

Conclusions

In this investigation we evaluate and quantify the ability of the SP model to predict the water depths, flooded area, and the hazard index using different spatial resolutions. The analysis is critical on the definition of the spatial scale required to develop reliable FRA and fast response early-warning systems.

The results yield important insights on the adequate resolutions required for the use of the PNSWE over NSWE considering the large reduction in computational time, analyzing the reproduced spatial distributions, and having a quantification of the errors induced by mesh coarsening.

Here we developed a SP model computational implementation in a Cartesian frame, introducing the building topography information at a sub-grid scale to reduce the computational time. To test the model with the most severe flooding conditions, we simulated dam-break flows for an idealized non-periodic city and a real-world urban zone in La Reina district using different spatial resolutions. With the outputs of water depth and velocity, we assessed the ability to capture the maximum water depth, flood extension, and the hazard index according to the definition of [START_REF] Garcia | Debris flows of december 1999 in venezuela[END_REF].

We compared the SP results with a fine-scale reference model in order to estimate the error of the sub-grid scale parameterization.

In general, the SP model captures all the main features of the flow when the resolution is equal to or below the street characteristic width SW. This is observed in both the idealized and real-world case. Moreover, when the resolution is very close to SW, the reduction in computation time is considerable, taking about 5% of the classical model CPU time while the main features of the flow are captured. Thus, the model has been found to be of great value when used at the meso-scale for these resolutions. The time reductions highlight the significant potential of the SP model for real-time forecasting and to implement early-warning systems.

When using SP, however, we found a diffusive-like behavior, especially when increasing significantly the mesh resolution. This effect is mainly due to the fact that by coarsening the mesh, the resolution of obstacles through the porosity term is gradually lost, and the flow has fewer impediments to spread throughout the floodplain. As a consequence, the maximum water depth within the streets is generally underestimated by SP. Even though the diffusive-like behavior suggests that the flooded area will always be larger with SP, this result was not observed in the case 2 of La Reina district. The SP simulations obtained for this case showed a total area that is underestimated, as the grid cannot resolve details of the terrain that control the flood propagation. The area is underestimated since preferential flowpaths that are characteristic of the real urban area are not captured for coarse resolutions, many areas are not flooded for coarser meshes, and thus they present no hazard.

Results show that even though variables such as velocity may lose physical meaning when using a coarse mesh, composed variables for hazard assessment can maintain the main categories of the flood for resolutions below SW.

It is important to note that there is a tendency for SPOR to follow the bottom slope as the resolution increases too much compared with SW (like the 64m mesh for La Reina). This is explained because the porosity tends to be more homogeneous and the gradient is not able to capture the preferential flowpaths.

Finally, we obtained the smallest error in the idealized case compared to the real-world case, even at the same spatial resolution. This suggests that the error is directly affected by the complexity of the urban patterns and that it can potentially enhance the underestimation of some variables, like flow depth.

This emphasizes the need for testing these type of sub-grid models in realworld urban environments, where the meso-scale can deliver fairly accurate results with a much faster computation time than classical models, carefully help to develop fast-response models approaching real-time forecasting for proper flood management and efficient decision-making.

  aligned in the same direction and perpendicular to the main flow, which is used to test the performance of the model under different grid sizes without the interference of other factors. In the second case, we apply the model to real conditions in a densely populated urban area, which shows the ability of the model to capture the flood propagation at different scales in a complex urban geometry. considered, is used to evaluate the performance of the model in a simplified distribution. The urban zone, presented in Figure 1(a), incorporates diverse building shapes and sizes, a typical Street Width (SW) of 9m, and voids representing parks or open areas. A constant slope of -0.573 • is also imposed in the direction of the flow.

Figure 1 :

 1 Figure 1: Idealized urban domain, with orthogonal street distribution, aligned with the main flow. The images show the maximum resolution of 1m (a), and the downgraded spatial grid with a size of 12m (b).

Figure 2 :

 2 Figure 2: Quebrada de Ramón watershed in red. The yellow shaded marks the total inundated area of the 1993 event, and in green the section considered in this investigation (Data from Google Earth (c)).

Figure 3 :

 3 Figure 3: Computational domain for flood simulations in La Reina district. (a) Street distribution in the real urban area at a resolution of 1m; (b) Downgraded domain at a resolution of 32m.

Figure 4 :

 4 Figure 4: Comparison of the Mean Absolute Error (MAE) and Mean Bias Error (MBE) for the single porosity (SP) model in the two studied cases for different spatial resolutions.

Figure 6 :

 6 Figure 6: Flooded area for case 2 for different spatial resolutions. Grid sizes larger than the street width fail to reproduce the inundated area.

Figure 7 :

 7 Figure 7: Flooded area for case 1 for different resolutions shows that at least for 6 times the size of the original mesh, the SP models captures with good agreement the inundated region.

Figure 8 :

 8 Figure 8: Hazard regions for case 1 according to the classification by[START_REF] Garcia | Debris flows of december 1999 in venezuela[END_REF]. The systematic coarsening of the grid shows that most of the high-resolution characteristics are preserved up to the resolution of 6m.

Figure 9 :

 9 Figure 9: Hazard regions for case 2 according to the classification by Garcia and Lopez

Table 1 :

 1 Computational times and averaged time-step for case 1 at different spatial resolutions.

	Model Resolution [m] CPU time [s] %CPU time Mean ∆t [s]
	NSWE	1	1246	100%	0.2
	SP	2	152	12.4%	0.27
	SP	3	55	4.2%	0.39
	SP	4	25.5	2.2%	0.5
	SP	6	9.6	0.7%	0.66
	SP	12	1.7	0.1%	1.54

Table 2 :

 2 Computational times and averaged time-step for case 2 at different spatial resolutions.

	Model Resolution [m] CPU time [s] %CPU time Mean ∆t [s]
	NSWE	1	130911	100%	0.002
	SP	4	21294	16.27%	0.18
	SP	8	6580	5.03%	0.172
	SP	16	1325	1.01%	0.122
	SP	32	97	0.07%	0.410
	SP	64	15	0.04%	0.105

Acknowledgements

This work was supported by Inria through the FLOTTE associated team, funded by the international relations department. Additional support has been provided by Fondecyt grant 1191785, the supercomputing infrastructure of the NLHPC (ECM-02) and ANID/Fondap 2022 Grant 1522A0005

suggesting again that the SP model has a more diffusive behavior.

In Figure 6 we show the flood map in the real urban area with an arbitrary orientation of the streets. The simulation evidences problems with the flow directionality induced by the SP model. In the fine-scale simulation in Figure 6a), the north section of the urban zone became flooded, despite the topography showing steeper slopes in other directions. This finding strongly suggests that the flooding in that area was influenced by the presence of preferential flow paths induced by the orientation of the streets. In contrast, for the coarser mesh of 64m in Figure 6f), the flooding pattern follows the


domain: math stat

text:
Hans-Christian Herbig

Daniel Herden
email: danielherden@baylor.edu

Christopher Seaton
email: seatonc@rhodes.edu

A SYMMETRIC FUNCTION APPROACH TO POLYNOMIAL REGRESSION

Keywords: 2020 Mathematics Subject Classification. Primary 62J02; Secondary 05E05, 65F05 polynomial regression, Schur polynomials, Vandermonde determinants. C

We give an explicit solution formula for the polynomial regression problem in terms of Schur polynomials and Vandermonde determinants. We thereby generalize the work of Chang, Deng and Floater to the case of model functions of the form n i=1 a i x d i for some integer exponents d 1 > d 2 > . . . > dn ≥ 0 and phrase the results using Schur polynomials. Even though the solution circumvents the well-known problems with the forward stability of the normal equations it is only of practical value if n is small since the number of terms in the formula grows rapidly with the number m of data points. The formula can be evaluated essentially without rounding.

Introduction

Linear regression goes back to ideas of A.M. Legendre [START_REF] Legendre | Nouvelles méthodes pour la détermination des orbites des comètes[END_REF] and C.F. Gauß [START_REF] Friedrich | Theoria combinationis observationum erroribus minimis obnoxiae[END_REF], while the first work on nonlinear regression is from 1815 and is due to J.D. Gergonne [START_REF] Diez | The application of the method of least squares to the interpolation of sequences[END_REF]. Nowadays regression is implemented in most of the mathematical software as a numerical routine and is widely used in statistics and science. The basic idea is to fit a model function, depending linearly on some parameters, to a set of data points in such a way that the sum of squares of the errors of the approximation is minimal.

In this paper we discuss the univariate case where the model function is a polynomial of the following form. Let d = (d 1 , d 2 , . . . , d n ) be an n-tuple of integers such that d 1 > d 2 > . . . > d n ≥ 0. By a polynomial of type d we mean a univariate function of the form f (x) = n i=1 a i x di with a i ∈ C. Let us assume that we have data points x = (x 1 , x 2 , . . . , x m ) T and y = (y 1 , y 2 , . . . , y m ) T ∈ C m . The interpolation problem f (x i ) = y i for i = 1, . . . , m is equivalent to the linear system

     x d1 1 x d2 1 . . . x dn 1 x d1 2 x d2 2 . . . x dn 2 . . . . . . . . . x d1 m x d2 m . . . x dn m           a 1 a 2 . . . a n      =      y 1 y 2 . . . y m      . (1.1)
Usually, this system is overdetermined (m > n) and inconsistent. Let us denote the coefficient matrix of this system by A ∈ C m×n and use the shorthand notation Aa = y for (1.1). Then the associated normal equation system

A * Aa = A * y, (1.2)
is consistent, where A * = A T denotes the Hermitian transpose of A. In fact, a = (a 1 , a 2 , . . . , a n ) T is a solution to (1.2) if and only if f (x) = n i=1 a i x di minimizes the distance m i=1 |f (x i )y i | 2 = ||Aa -y|| 2 = d(Aa, y) 2 . Here we use the Hermitian inner product on C n . For d = (n -1, n -2, . . . , 1, 0) this is known as the problem of polynomial regression of degree n -1. The resulting minimal distance d(Aa, y) is given by ||y|| 2 -A * y|a . Note that the coefficient matrix of (1.2),

A * A = m l=1 x di l x dj l i,j=1,...,n , (1.3)
is the Hankel matrix of the sequence of power sums if x ∈ R m and d = (n -1, n -2, . . . , 1, 0).

In this paper we exploit the fact that the system (1.2) is equivariant with respect to simultaneous permutation of the x i s and the y i s. It has a unique solution if and only if the columns of A are linearly independent, and we will use symmetric polynomials to provide concrete expressions (see Theorem 2.2 below) for this unique solution of (1.2). If d = (n -1, n -2, . . . , 1, 0) the solution is unique if and only of the data points are distinct. In fact we will construct for an injective A a matrix B ∈ C n×( m n-1 ) such that the solution operator A + to the minimal squares problem (also known as the Moore-Penrose pseudoinverse) can be written as

A + = BB * A * = B(AB) * . (1.4)
The projection to the image of A is then P = AB(AB) * . Writing the complementary projection as P ⊥ = 1 -P = 1 -AB(AB) * the error of the regression problem can be expressed as d(Av, u) = u|P ⊥ u . The formula for d = (n -1, n -2, . . . , 1, 0) was already elaborated previously by Chang, Deng and Floater in [START_REF] Chang | An interpolatory view of polynomial least squares approximation[END_REF] for weighted least squares without making use of symmetric functions. The weighted case will be discussed in Section 3.

At this stage it is already apparent that the approach has practical limitations. The problem is that the number of columns of B is From the point of view of stability, however, the situation is much better. The least squares problem has a forward error that is measured by κ 2 (A) 2 if the error of the minimal squares approximation is large and κ 2 (A) if it is small (see [START_REF] Nicholas | Accuracy and stability of numerical algorithms[END_REF]Theorem 20.1]). Here κ 2 (A) denotes the condition number of the matrix A1 . The forward error for solving the normal equation directly is, on the other hand, about κ 2 (A) 2 (see [START_REF] Nicholas | Accuracy and stability of numerical algorithms[END_REF]Subsection 20.4]). Therefore, from the point of view of forward stability, the normal equation are considered problematic. This is particularly important for regression since κ 2 (A) for matrices of Vandermonde type can be large (see, for example, [START_REF] Victor | How bad are Vandermonde matrices?[END_REF]). Our explicit solution does not have this defect. Problems with the stability can only occur in the polynomial evaluations and matrix multiplications involved in the formula. But in principle those calculations can be done essentially without rounding (only at the last step there is a division).

The solution formula

We now recall the definition of Schur polynomials. Any decreasing sequence λ = (λ 1 , . . . , λ n ) of integers λ 1 ≥ λ 2 ≥ . . . ≥ λ n ≥ 0 represents a partition of N = n i=1 λ i into at most n parts. We denote by δ(n) the partition (n-1, n-2, . . . , 1, 0) of N = n 2 . We introduce the alternating polynomial

a λ (z) = |(z λj i ) i,j=1,...,n | for z = (z 1 , . . . , z n ). Note that a δ(n) (z) = 1≤i<j≤n (z i -z j ) =: V (z) is the Vandermonde determinant. The Schur polynomial associated to λ = (λ 1 , . . . , λ n ) is s λ (z) := a λ+δ(n) (z) a δ(n) (z) = |(z λj +n-j i ) i,j=1,...,n | V (z)
,

where the partition λ + δ(n) results from the componentwise addition of λ and δ(n). This, indeed, defines a polynomial as every alternating polynomial is divisible by the Vandermonde determinant. In fact, s λ is a polynomial with non-negative integer coefficients that can be determined by counting out semistandard Young tableaux of shape λ.

Given x := (x 1 , . . . , x m ) ∈ C m and k = (k 1 , k 2 , . . . , k n ) with 1 ≤ k 1 < k 2 < . . . < k n ≤ m, we write x k := (x k1 , . . . , x kn ) ∈ C n .
We can view k as an element of [m] n , the collection of all n-element subsets of [m] = {1, 2, . . . , m}. We note the following immediate result. Proposition 2.1. For all m, n ∈ Z >0 , and partitions λ = (λ 1 , . . . , λ n ) and ν = (ν 1 , . . . , ν n ), and

x 1 , x 2 , . . . , x m ∈ C holds m k=1 x λi+n-i k x νj +n-j k i,j=1,...,n = k∈( [m] n ) s λ (x k )s ν (x k )|V (x k )| 2 .
Proof. With the Cauchy-Binet formula we have

m k=1 x λi+n-i k x νj +n-j k i,j=1,...,n = (x λi+n-i k ) T k∈[m],i∈[n] (x νj +n-j k ) k∈[m],j∈[n] = 1≤k1<...<kn≤m (x λj +n-j ki ) T i,j=1,...,n (x νj +n-j ki ) i,j=1,...,n = 1≤k1<...<kn≤m s λ (x k )V (x k )s ν (x k )V (x k ) = k∈( [m] n ) s λ (x k )s ν (x k )|V (x k )| 2 .
We associate to d

= (d 1 , . . . , d n ) a partition λ = (λ 1 , . . . , λ n ) by λ = d -δ(n), i.e., λ k := d k + k -n for k = 1, . . . , n. Note that λ is a partition of D -n 2 . Moreover, let d i = (d 1 , . . . , d i-1 , d i+1 , . . . , d n ) and put λ[i] = d i -δ(n -1). We see that λ[i] is a partition of D -d i -n-1
2 . We apply Proposition 2.1 to obtain the following result.

Theorem 2.2. Let x 1 , . . . , x m ∈ C with m ≥ n such that the system (1.2) has a unique solution.

Then, using the notation above, this unique solution a = (a 1 , . . . , a n ) as a function of d, x, y is

a i = n j=1 (-1) i+j   l∈( [m]
n-1 )

s λ[i] (x l )s λ[j] (x l )|V (x l )| 2   m k=1 x dj k y k k∈( [m] n ) |s λ (x k )V (x k )| 2 (2.1)
for i = 1, . . . , n. In particular, if x 1 , . . . , x m are positive real numbers with at least n distinct values, then (1.2) has a unique solution.

Proof. The system (1.2) has a unique solution if and only if the matrix A * A is invertible, and we have a = (A * A) -1 A * y. Applying the formula for the inverse of A * A gives

a i = n j=1 (-1) i+j M ji (A * y) j |A * A| ,
where M ji denotes the determinant of the submatrix of A * A formed by deleting the jth row and ith column of

A * A. We have (A * y) j = m k=1
x dj k y k for the jth component of the vector A * y. Moreover, combining (1.3) and Proposition 2.1 gives

|A * A| = m k=1 x λi+n-i k x λj +n-j k i,j=1,...,n = k∈( [m] n ) s λ (x k )s λ (x k )|V (x k )| 2 = k∈( [m] n ) |s λ (x k )V (x k )| 2 .
A similar argument applies for

M ji = l∈( [m] n-1 ) s λ[i] (x l )s λ[j] (x l )|V (x l )| 2 .
Finally, recall that s λ is a polynomial with non-negative integer coefficients. Thus, V (x k ) is nonzero for any vector x k of distinct complex numbers, and s λ (x k ) is positive for any vector x k of positive real numbers.

We can reformulate Theorem 2.2 also in terms of the solution operator A + = BB * A * of the minimal squares problem, i.e., the Moore-Penrose pseudoinverse. The matrix B ∈ C n×( m n-1 ) is given by the formula,

B il := (-1) i a λ[i] (x l ) k∈( [m] n ) |a λ (x k )| 2 = (-1) i s λ[i] (x l )V (x l ) k∈( [m] n ) |s λ (x k )V (x k )| 2 , (2.2)
where we understand the set [m] n-1 to be endowed with some fixed total order for indexing the columns of the matrix B.

2.1. The case of polynomial regression. Consider the case d = (n-1, n-2, . . . , 1, 0) of polynomial regression of degree n -1. In this case, λ k = 0 for all k so that λ

[i] = (1, . . . , 1 i-1 , 0, . . . , 0 n-i ).
As s λ = 1, Equation (2.1) becomes

a i = n j=1 m k=1 (-1) i+j x n-j k y k l∈( [m] n-1 ) s λ[i] (x l )s λ[j] (x l )|V (x l )| 2 k∈( [m] n ) |V (x k )| 2 . (2.3)

Weighted regression

What has been said in Section 2 can be easily adapted to regression with weighted least squares (see [START_REF] Chang | An interpolatory view of polynomial least squares approximation[END_REF]). Let us first investigate what happens to the normal equations if we work with a non-standard inner product. Proof. Let G = W * W be the Gram matrix so that u|v := u * Gv. We have u * (A )

* x = u * GAx for all u ∈ C m , x ∈ C n so that (A ) * = GA.
This in turn can be written as A = (GA) * and hence

A Aa = A * G * Aa = (W A) * W A a and A * G * y = (W A) * W y.
We now assume that W = diag(w 1 , . . . , w m ) for some weights w 1 , . . . , w m ∈ C × . Instead of minimizing Aa -y, Aa -y we are now going to minimize Aa -y|Aa -y , i.e., to solve (W A) * W A a = (W A) * W y. We introduce the shorthand notation

|w k | 2 := n j=1 |w kj | 2 ,
where k ∈ [m] n is interpreted as 1 ≤ k 1 < k 2 < . . . < k n ≤ m. A little modification of the argument of Proposition 2.1 yields the formula

m k=1 |w k | 2 x λi+n-i k x νj +n-j k i,j=1,...,n = k∈( [m] n ) |w k | 2 s λ (x k )s ν (x k )|V (x k )| 2 .
As a consequence the solution formula (2.1) of Theorem 2.2 becomes

a i =   l∈( [m] n-1 ) |w l | 2 s λ[j] (x l )s λ[j] (x l )|V (x l )| 2   m k=1 |w k | 2 x dj k y k k∈( [m] n ) |w k | 2 s λ (x k )s ν (x k )|V (x k )| 2 ,
while the pseudoinverse turns into BB * A * W * W with

B il = (-1) i w l s λ[i] (x l )V (x l ) k∈( [m] n ) |w k | 2 |s λ (x k )V (x k )| 2 .
s λ (x 1 , x 2 , x 3 ) = (x 1 + x 2 )(x 1 + x 3 )(x 2 + x 3 ), s λ[1] (x 1 , x 2 ) = x 1 + x 2 , s λ[2] (x 1 , x 2 ) = (x 1 + x 2 )(x 2 1 + x 2 2 ), s λ[3] (x 1 , x 2 ) = x 2 1 x 2 2 (x 1 + x 2 )
. Then Equation (2.1) yields

a 1 = 3 j=1 (-1) j+1 l∈( [m]
2 )

s λ[1] (x l )s λ[j] (x l )|V (x l )| 2 m k=1 x dj k y k k∈( [m]
3 )

|s λ (x k )V (x k )| 2 = l∈( [m] 2 ) |x 2 l1 -x 2 l2 | 2 m k=1 x 4 k y k -(x 2 l1 + x 2 l2 )x 2 k y k + x 2 l1 x 2 l2 y k k∈( [m] 3 ) (x 2 k1 -x 2 k2 )(x 2 k1 -x 2 k3 )(x 2 k2 -x 2 k3 ) 2 , a 2 = 3 j=1 (-1) j l∈( [m]
2 )

s λ[2] (x l )s λ[j] (x l )|V (x l )| 2 m k=1 x dj k y k k∈( [m]
3 )

|s λ (x k )V (x k )| 2 = - l∈( [m]
2 )

|x 2 l1 -x 2 l2 | 2 (x 2 l1 + x 2 l2 ) m k=1 x 4 k y k -(x 2 l1 + x 2 l2 )x 2 k y k + x 2 l1 x 2 l2 y k k∈( [m] 3 ) (x 2 k1 -x 2 k2 )(x 2 k1 -x 2 k3 )(x 2 k2 -x 2 k3 ) 2 ,
and

a 3 = 3 j=1 (-1) j+1 l∈( [m]
2 )

s λ[3] (x l )s λ[j] (x l )|V (x l )| 2 m k=1 x dj k y k k∈( [m]
3 )

|s λ (x k )V (x k )| 2 = l∈( [m]
2 )

|x 2 l1 -x 2 l2 | 2 x 2 l1 x 2 l2 m k=1 x 4 k y k -(x 2 l1 + x 2 l2 )x 2 k y k + x 2 l1 x 2 l2 y k k∈( [m] 3 ) (x 2 k1 -x 2 k2 )(x 2 k1 -x 2 k3 )(x 2 k2 -x 2 k3 ) 2 .
The matrix entry B il is given by .

B il = (-1) i (x l1 -x l2 )s λ[i] (x l ) k∈( [m] 3 ) (x 2 k1 -x 2 k2 )(x 2 k1 -x 2 k3 )(x 2 k2 -x 2 k3 ) 2 , i.e., B 1,l = x 2 l2 -x 2 l1 k∈( [m] 3 ) (x 2 k1 -x 2 k2 )(x 2 k1 -x 2 k3 )(x 2 k2 -x 2 k3 ) 2 , B 2,l = x 4 l1 -x 4 l2 k∈( [m] 3 ) (x 2 k1 -x 2 k2 )(x 2 k1 -x 2 k3 )(x 2 k2 -x 2 k3 ) 2 , B 3,l = x 2 l1 x 2 l2 (x 2 l2 -x 2 l1 ) k∈( [m] 3 ) (x 2 k1 -x 2 k2 )(x 2 k1 -x 2 k3 )(x 2 k2 -x 2

Complexity

Here we are examining the complexity of possible evaluations of Equation (1.4). The simple message is that, for fixed n and d 1 , it can be estimated asymptotically as O (m n ) for m → ∞ (for an introduction to asymptotic analysis the reader may consult [START_REF] Graham | Concrete mathematics: a foundation for computer science[END_REF]). The dominating contributions are coming from the matrix multiplication and the evaluation of the denominator. The coefficients of the estimates depend on n and d 1 , the details of these dependencies depend on the details of the evaluation. The presentation the material is hopefully intelligible for people that are not accustomed to counting floating point operations (FLOPs). We have to conclude that for n > 3 evaluating our formulas is not competitive with the numerical algorithms commonly used.

First, recall that the complexity of naively multiplying three matrices

X ∈ C n×m , Y ∈ C m×k , Z ∈ C k×l is O(nmk + mkl). So the complexity of naively multiplying BB * A * is O n(n + m) m n-1 as B ∈ C n×( m n-1 ) , B * ∈ C ( m n-1 )×n , A * ∈ C n×m . Since for positive integers M, N we have M N N ≤ M N ≤ eM N N we can estimate n(n + m) m n -1 ≤ n(n + m) e n -1 n-1 m n-1 .
m n ∼ n 2 e em n
n . This is because

n n = (n -1) n n n -1 n = (n -1)(n -1) n-1 1 1 -1/n n ∼ (n -1)(n -1) n-1 e ,
which implies

1 (n-1) n-1 ∼ (n-1) en n . From this, we deduce O n(n + m) m n-1 = O n 2 e em n
n . Now we would like to investigate the complexity of evaluating the matrix B of Equation (2.2). In the denominator we have to evaluate s λ (x k ). We could do that using a multivariate Horner's scheme for evaluating a polynomial f , the number of FLOPs being 2

deg(f )+m m -2 (see [2, Theorem 3.1]). Hence with deg(s λ (x k )) = |λ| = n k=1 (d k +k -n) = n+1 2 -n 2 + n k=1 d k = -n 2 + n k=1 d k this is 2 |λ|+m m -2. But 2 |λ|+m m -2 = O ((|λ| + m) m
) grows extremely rapidly with the number of data points m. A better choice seems to be the second Jacobi-Trudi identity (see [START_REF] Macdonald | Symmetric functions and Hall polynomials[END_REF]Equation (3.5)]), expressing s λ (x k ) as a λ 1 × λ 1 -determinant in the elementary symmetric polynomials. To evaluate the elementary symmetric polynomials one has to expand (X -

x 1 )(X -x 2 ) . . . (X -x m ). If (X -x 1 )(X -x 2 ) . . . (X -x k-1 ) is already expanded, then expanding (X -x 1 )(X -x 2 ) . . . (X -x k-1 ) • (X -x k ) just needs 2(k -1)
extra multiplications and k extra additions, i.e., 3(k -1) + 1 more FLOPs. Hence evaluating the elementary symmetric polynomials costs

m k=1 (3(k -1) + 1) = 3 m 2 + m = O m 2
FLOPs. Calculating the determinant costs in the worst case O λ 3 1 FLOPs (see [START_REF] Nicholas | Accuracy and stability of numerical algorithms[END_REF]). For the specific Schur polynomials in question λ 1 = d 1 , and the cost of evaluating such a Schur polynomial can be estimated as O m 2 + d

O m 2 + m n d 3 1 = O m 2 + e n n m n d 3 1 .
For m n ≥ 2 the cost of it is roughly the same as that of the matrix multiplication. The numerator can be analyzed along similar lines, yielding an estimate

O m 2 + n m n -1 d 3 1 = O m 2 + n e n -1 n-1 m n-1 d 3 1 .
Note here that the evaluations of the elementary symmetric polynomials can be repurposed, hence they contribute m(3m + 2) = O(m 2 ) FLOPs. For m n > 2 the cost of evaluating the numerator is negligible in comparison to the contribution of the denominator. 5.1. A numerical example. Let us look into a concrete example of the evaluation that we implemented into Mathematica [10] for d = (4, 2, 0), i.e., the model function to fit is again

f (x) = a 1 x 4 + a 2 x 2 + a 3 .
To give a practical application, this could, for example, be taken as a model function used for digital image processing of a human mouth. Let us assume for the moment f (x) to be non-constant. From f (x) = 4a 1 x 3 + 2a 2 x = 2x 2a 1 x 2 + a 2 we deduce that for a 1 a 2 ≥ 0 there is just one extremum. Otherwise there are three extrema x = 0 and x = ± -a 2 /(2a 1 ). Here we take the liberty to ignore the non-generic case a 1 = 0. From the sign of a 1 a 2 a machine could infer if the mouth is smiling (i.e., a 1 a 2 > 0) or frowning (i.e., a 1 a 2 < 0) 2 .

To investigate empirically the running time of the evaluation of our Formula (2.1) we created data points by superimposing random noise on the evaluations of The range for x is here [-500, 500]∩ Z. In Figure 1 we show an example of a regression obtained this way with = 10, i.e., m = 101 data points. In Figure 2 we depict how the running time of our implementation rises when we increase the number of data points.

Recursive evaluations

Observe that, when fitting a fixed polynomial f (x) = n i=1 a i x di to a large data set or multiple data sets, the partitions λ and λ[i], i = 1, . . . , n, depend only on d. Therefore, the polynomials s λ[i] and s λ can be computed and stored once (e.g., using the second Jacobi-Trudi identity). In addition, if a data set is enlarged, Equation (2.1) can be used recursively as we now describe. In this way, amending one additional data point to a data set of m -1 data points reduces the complexity of the evaluations from O(m n ) to O(m n-1 ) for large m.

Suppose that one has applied Equation (2.1) to a collection of data points x = (x 1 , x 2 , . . . , x m ) T and y = (y 1 , y 2 , . . . , y m ) T . In this process, one stores not only the a i , but the following quantities. For i = 1, . . . , n and j = 1, . . . , n, define

S i,j (x) = l∈( [m] n-1 ) s λ[i] (x l )s λ[j] (x l )|V (x l )| 2 , T j (x, y) = m k=1 x dj k y k ,
and let

N i (x, y) = n j=1 (-1) i+j S i,j (x)T j (x, y) and D(x) = k∈( [m] n ) |s λ (x k )V (x k )| 2
denote the numerators and denominator in the expression for a i in Equation (2.1). Enlarging the data set by one point, let x = (x 1 , x 2 , . . . , x m , x m+1 ) T and y = (y 1 , y 2 , . . . , y m , y m+1 ) T We express N i (x , y ) as

n j=1 (-1) i+j    l∈( [m+1] n-1 ) s λ[i] (x l )s λ[j] (x l )|V (x l )| 2    m+1 k=1 x dj k y k = n j=1 (-1) i+j   Si,j(x, y) + l∈( [m] n-2 ) s λ[i] (x l,m+1 )s λ[j] (x l,m+1 )|V (x l,m+1 )| 2    × T j (x, y) + x dj m+1 y m+1 = N i (x, y) + n j=1 (-1) i+j S i,j (x, y)x dj m+1 y m+1 + l∈( [m] n-2 ) s λ[i] (x l,m+1 )s λ[j] (x l,m+1 )|V (x l,m+1 )| 2 l∈( [m] n-2 ) T j (x, y) + x dj m+1 y m+1 . Let R i,j (x , y ) = l∈( [m] n-2 ) s λ[i] (x l,m+1 )s λ[j] (x l,m+1 )|V (x l,m+1 )| 2 ,
and then we can express N i (x , y ) as

N i (x, y) + n j=1 (-1) i+j S i,j (x, y) + R i,j (x , y ) x dj m+1 y m+1 + R i,j (x , y )T j (x, y) .
Similarly, D(x, x m+1 ) can be written

k∈( [m+1] n ) |s λ (x k )V (x k )| 2 = D(x) + k∈( [m] n-1 ) |s λ (x k,m+1 )V (x k,m+1 )| 2 .
This yields the following expression for the a i to fit to the data set (x , y ):

(6.1) a i = Dai + n j=1 (-1) i+j (Si,j + Ri,j) x d j m+1 ym+1 + Ri,jTj D + k∈( [m] n-1 ) |s λ (x k,m+1 )V (x k,m+1 )| 2
for i = 1, . . . , n, where D = D(x), S i,j = S i,j (x, y), etc., and R i,j = R i,j (x , y ).

6.1. Recursive evaluation of even polynomials of degree 4. As an illustration, let us consider the case of fitting a polynomial of the form f (x) = a 1 x 4 + a 2 x 2 + a 3 treated in Section 4.1.

Fitting such a polynomial to a data set x = (x 1 , x 2 , . . . , x m ) T and y = (y 1 , y 2 , . . . , y m ) T involves computing

S 1,1 (x) = l∈( [m]
2 )

|x 2 l1 -x 2 l2 | 2 , S 1,2 (x) = S 2,1 (x) = l∈( [m]
2 )

(x 2 l1 + x 2 l2 )|x 2 l1 -x 2 l2 | 2 , S 1,3 (x) = S 3,1 (x) = l∈( [m]
2 ) 2 ) 3 )

x 2 l1 x 2 l2 |x 2 l1 -x 2 l2 | 2 , S 2
x 2 l1 x 2 l2 (x 2 l1 + x 2 l2 )|x 2 l1 -x 2 l2 | 2 , S 3
(x 2 k1 -x 2 k2 )(x 2 k1 -x 2 k3 )(x 2 k2 -x 2 k3 ) 2 .
Given these quantities, one can fit the data set x = (x, x m+1 ) and y = (y, y m+1 ) with an additional point by computing 2 )

(x 2 k1 -x 2 k2 )(x 2 k1 -x 2 m+1 )(x 2 k2 -x 2 m+1 ) 2
as well as 3m additional entries B i,l,m+1 for i = 1, 2, 3 and l = 1, . . . , m.

m n- 1 ,

 1 which is only reasonable in the typical situation when m n for roughly the values n = 2, 3, 4, 5. The number of floating point operations (FLOPs) for the naive multiplication BB * A * is then n(n + m) m n-1 . For m n 0 this number of operations can be estimated by n(n + m) m n-1 = O n 2 e em n n . Similar estimates appear when counting the floating point operations for evaluating B (see Section 5).

Lemma 3 . 1 .

 31 Let W ∈ GL m (C) be an invertible matrix and define the inner product u|v := u * W * W v on C m . Let A denote the adjoint of A ∈ C m×n = Hom (C n , C m ) with respect to | , i.e., A u|x = u|Ax for all u ∈ C m , x ∈ C n . Then the normal equation A Aa = A y for a, y ∈ C n can be written as (W A) * W A a = (W A) * W y.

4 . Examples 4 . 1 .

 441 Even polynomials of degree 4. Consider the case d = (4, 2, 0) of fitting an even polynomial of the form f (x) = a 1 x 4 + a 2 x 2 + a 3 . In this case, λ = (2, 1, 0); λ[1] = (1, 0); λ[2] = (3, 0); and λ[3] = (3, 2); from which one computes

k3 ) 2 . 4 . 2 .

 242 Fitting power functions. Suppose d = (d), i.e., the function to fit is the power function f (x) = x d . Then λ = λ = d so that s λ = x d , and d 1 = λ[START_REF] Chang | An interpolatory view of polynomial least squares approximation[END_REF] is the empty partition. Equation (

For n m this in turn can be estimated by O nm m n- 1 =

 1

3 1 .

 1 Evaluating the full Vandermonde determinant costs m 2 = O(m 2 ) FLOPs, and hence evaluating all of those in the denominator by diving out linear factors m 2 + 2m(m -n) = O(m 2 ) FLOPs. So the complexity of evaluating the denominator is

2

  In physics this type of model functions plays a role in the study of phase transitions associated to the phenomenon of spontaneous symmetry breaking.

Figure 1 .

 1 Figure 1. Here = 10. The dashed graph is the one the function g(x) of Equation (5.1), the solid one is the graph of f (x).

Figure 2 .

 2 Figure 2. Here we put a doubly logarithmic plot of the running time (measured in seconds) of our implementation on the y-axis against the number of data points on the x-axis. The slope equals approximately 3, verifying empirically our predicted complexity of the evaluation as O(m 3 ).

-x 4 l2 | 2 ,S 2 , 3

 4223 (x) = S 3,2 (x) = l∈( [m]

2 )

 2 ,3 (x) = l∈( [m] |x l1 | 4 |x l2 | 4 |x 2 l1 -x 2 l2 | 2 , T j (x, y) = m k=1 x dj k y k , j = 1, 2, 3; (d 1 , d 2 , d 3 ) = (4, 2, 0), N i (x, y) = (-1) i -S i,1 (x)T 1 (x, y) + S i,2 (x)T 2 (x, y) -S i,3 (x)T 3 (x, y) , i = 1, 2, 3,andD(x) = k∈( [m]

R 1 , 1 (x , y ) = m l=1 |x 2 l -x 2 m+1 | 2 ,R 1 , 2 (x 2 l + x 2 m+1 )|x 2 l -x 2 m+1 | 2 ,R 1 , 3 l=1 x 2 l x 2 m+1 |x 2 l -x 2 m+1 | 2 ,R 2 , 2 (x , y ) = m l=1 |x 4 l -x 4 m+1 | 2 , 2 l2 -x 2 l1D 3 .

 112221222221322222442223 (x , y ) = R 2,1 (x , y ) = m l=1 (x , y ) = R 3,1 (x , y ) = m R 2,3 (x , y ) = R 3,2 (x , y ) = m l=1 x 2 l x 2 m+1 (x 2 l + x 2 m+1 )|x 2 l -x 2 m+1 | 2 , R 3,3 (x , y ) = m l=1 |x l | 4 |x m+1 | 4 |x 2 l -x 2 m+1 | 2 ,as well as the new denominator terms k∈( the data set (x, y), one may remember the matrix B as the denominator D(x) above and the entriesB 1,l = xTo compute the new matrix B for fitting an enlarged data set, one would compute the new denominator D(x) = k∈([m]

It can be expressed, for example, as the ratio of the biggest over the smallest singular value of A.

Acknowledgments C.S. was supported by an AMS-Simons Research Enhancement Grant for PUI Faculty and a Rhodes College Faculty Development Grant.

for PUI Faculty and a Rhodes College Faculty Development Grant.