Amine Ammar
email: amine.ammar@ujf-grenoble.fr

Francisco Chinesta
email: francisco.chinesta@paris.ensam.fr

David Ryckelynck
email: david.ryckelynck@paris.ensam.fr

Deterministic particle approach of Multi Bead-Spring polymer models

Keywords: Fokker-Planck equation, Multi Bead Spring models, method of particles. MOTS-CLÉS : équation de Fokker-Planck, modèles Multi-Haltères, méthodes particulaires

de documents 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.

Introduction

As a consequence of the increasing use of composite materials, micro-macro simulation comes back into fashion nowadays with the development of computer performance. However, it remains crucial that simplicity and performance of numerical techniques affect the usability of microscopic models.

In the case of dilute polymer, simple dumbbell models, such as the FENE model, are not representative of the behaviour of a polymeric chain. In fact, both, the flexibility of the chain and the description of its complex configurations state are better accounted using a Multi Bead-Spring (MBS) representation.

Simulation of MBS has been widely performed by using stochastic techniques, whose accuracy is related to the number of realizations considered. In the context of deterministic techniques, the use of standard discretization techniques such as finite differences or finite elements are forbidden, because of the extremely large number of degrees of freedom involved due to the multidimensional character of the associated Fokker-Planck equation. Thus, the direct resolution of the Fokker-Planck equation in the finite element framework only concerned simple dumbbell models [START_REF] Lozinski | A fast solver for Fokker-Planck equation applied to viscoelastic flows calculations: 2D FENE model[END_REF][START_REF] Chauviere | Simulation of dilute polymer solutions using a Fokker-Planck equation[END_REF]). An alternative to the standard discretization techniques of the Fokker-Planck equation consists of the use of the method of particles. However, this technique has never been used until now for solving models defined in a multidimensional space, as in the case of the MBS polymer models. The main aim of this work is the analysis of the application of this technique to those models.

Results of such simulation highlight the flexibility effects of the chain (especially in shearing flows), which cannot be observed in the case of simple dumbbell models. However, a number of difficulties persist, some of them restricting seriously its application in some simulation conditions.

Mechanical model

As shown in Figure 1, the MBS chain consists of N+1 beads connected by N springs. The bead serves as an interaction point with the solvent and the spring contains the local stiffness information depending on local stretching (see Bird et al., 1987, for more details).

Figure 1. MBS model for polymer chain

The dynamics of the chain is governed by hydrostatic, Brownian and connector forces. If we denote by j r the velocity of the bead 'j' and by j q the velocity of the spring connector, then we obtain:

N j r r q j j j ,.., 1 , 1 [1]
The dynamics of each bead can be written as: 01

Hydrostatic effects

Interactions Forces Brownian effects (g r a d ) l n ( ) 0

cc jj b j j j rv v r k T F F r [2]
where is the drag coefficient, v is the velocity field, v 0 is an average velocity, k b is the Boltzman constant, T is the absolute temperature and is the probability distribution function. From Equations [1] and [2] it results:

1 1 grad ln( ) N c jj j k b k k k qv q A k T F q [3]
where A jk is the Rouse Matrix defined by:

1 if 1 if 2 k j k j A jk
In the Rouse model connector force F c is a linear function of the connector orientation vector. But in our case a Finitely Extensible Non Linear spring will be used so that each spring branch behaves as in FENE model whose connector force will be detailed later.

Using now the evolution of the distribution function given by N j N j j N t q q q q t t q q 1 1 1 ) , ,.., ( ) , ,.., (

1 grad ( ) 4 c qv q A F q [6] 1 ( ) 4 qA tq q q =7]
Where q denotes the vector containing the N spring connector vectors j q , with 1 jN ,, .

The problem [6]- [7] is defined in a MxN-dimensional space, being M=1,2 or 3 (the dimension of the physical space) and N the number of springs that compose the polymeric chain. In the present work the connector force is given by:

j j j c q b q q F / 1 1 ) ( 2 [8]
where b is the maximum stretching of each spring connector of the chain.

Deterministic discretization approach using the method of particles

Advection form of the Fokker-Planck equation

For applying the method of particles to discretize the Fokker-Planck equation, we firstly introduce its diffusion term into the advection one (Ammar and Chinesta 2005;[START_REF] Chinesta | On the solution of the Fokker-Planck equations in steady recirculating flows involving short fiber suspensions[END_REF][START_REF] Chaubal | Smoothed particle hydrodynamics techniques for the solution of kinetic theory problems[END_REF]. Thus, it results the pseudoadvection field:

~1 4 q qq A [9]
being the associated Fokker-Planck equation:

( ) q tq [10] Problem [9]-[10] is equivalent to the initial problem [6]-[7].

Smoothed particle approximation

In order to discretize Equations [9]-[10] one must define the approximation of the unknown field which is, in our case, the distribution function. The method of particles uses a Dirac representation, consisting of ()

Cn i i i
qq , where the superscript refers to each chain of the population involved in the simulation that consists of Cn chains (particles). However, from a computational point of view the derivative of Dirac masses cannot be defined in a proper way. So, usually, a smoothed approximation of the Dirac mass is introduced according to: 2 ( ) with d d de d [11] where denotes a smoothing parameter and 'd' denotes a distance in the MxNdimensional space of molecule configurations.

Thus we can write:

( ) i d i i qq d [12]
where i i dq q is the Euclidean norm in the MxN-dimensional space.

The smoothed distribution can then be written as: 1 ()

Cn i i i qd [13]
Where i is the weight related to each chain, which, in this work, is considered the same overall the population ( = 1/Cn for satisfying the normality condition of the distribution function).

Equilibrium conditions

Before starting simulations, an equilibrium state eq can be defined, and then used for defining the initial condition for start-up problems, as well as for adjusting the value of the smoothing parameter by comparing eq with the solution computed from the steady state related to the transient solution of the Fokker-Planck equation without the advection term.

According to Bird et al., (1987) and after the application of the normality condition in the 2xN dimensional case considered in this work, the expression of eq (for the connector force defined in Equation [8] is given by:

N N j b j eq b b b q q 2 2 1 ) ( 1 2 / 2
[14]

Simulation algorithm

The simulation algorithm proceeds as follows: -Consider Cn particles (molecules) consisting of N connectors whose orientation must describe the initial known distribution: (0) ,1 ,,; 1 , ,

i j qi C n j N that is, ( 0
)
i q is assumed known.

Index i refers to the chain (particle) involved in the simulation, j to the connector in the chain, and (n) indicates the time step.

-For each time step: -Compute the advection field according to Equation [9]:

() () () ~() () () 1 4 () in n in q in in n q qq A q where () () 1 () Cn ni n i i qd , being () () in in dq q .
Update the molecular conformation:

() (1 ) ( ) in in in qqq t
Update the distribution function:

(1 ) (1 ) 1 ()

Cn ni n i i qd

Numerical test

In order to validate the approach just proposed we consider firstly a problem whose solution can be obtained analytically. This is for example the case of a single FENE dumbbell in a one-dimensional physical space. This problem was deeply analysed in [START_REF] Keunings | On the peterlin approximation for finitely extensible dumbbells[END_REF]. The number of particles involved in the simulation is in the order of 2400, being the other simulation parameters the following: 0.001 t [15]

that is equivalent to Equation [14] after normalization. We prefer using in this numerical example Equation [15] as initial condition instead of Equation [14] because in this manner the computed solution can be directly compared with the results reported in [START_REF] Keunings | On the peterlin approximation for finitely extensible dumbbells[END_REF].

The evolution of the distribution function is depicted in Figure 2 which represents the solutions at times: t=0,0.05,0.1,0.15,0.2,0.3,0.4,0.5,1.

Figure 2. Evolution of the distribution function related to a single 1D FENE model

Figure 3 zooms the region located around the maximum molecule stretching. These distributions are in perfect agreement with the exact solution as well as with the numerical solutions computed using the stochastic technique that can be found in [START_REF] Keunings | On the peterlin approximation for finitely extensible dumbbells[END_REF]. Figure 4 depicts the histogram related to the particle position at t=1 (we recall that the portion of the particle indicates the position of one of its ends, the other one remains at the origin x=0). Simular results were obtained by using only around 600 particles. Now, after proving the accuracy of the proposed strategy, we consider higher dimensional problems. For this purpose we consider a MBS -FENE model defined in a one-dimensional physical space composed of two dumbbells. In this case a reference solution can be easily obtained by using a more experienced numerical technique as is the finite element method. We consider the problem defined by the following parameters: Cn=500 (number of particles), b=10 ( b being the dimensionless parameter giving the maximum stretching), N = 2 (number of branches), M=1 (dimension of the physical space) and 0.35 du We dx

.

Figure 5 depicts the reference steady solution as well as the solution computed by using the method of particles with a smoothing parameter 0.4 . The steady distribution function depends on the stretching of both branches 1 q and 2 q .

Figure 5. Distribution function related to a 2 branches FENE MBS molecule model conformation: (left) Reference solution; (right) Solution computed using the method of particles

The distribution depicted in Figure 5(right) is related to the particle distribution depicted in Figure 6. In order to conclude about the incidence of the smoothing parameter on the computed solution we show in Figure 7 the distributions computed by using different values of that parameter, in particular 0.1,0.2,0.3,0.5 ( 0.4 being the optimum value whose associated distribution is depicted in Figure 5(right)). We can notice the significant impact of the smoothing parameter on the particles distribution and then on the conformation distribution. Moreover this parameter must be identified in each particular application because it depends on the number of particles considered in the simulation; on the solution itself (more localized is the solution smaller must be the smoothing parameter, but at the same time it must be large enough for coupling particles in the regions with low particles density); and on the dimension of the physical and conformation spaces. The definition of an appropriate adaptive procedure for evaluating the local value of the smoothing parameter is not an easy task. 0.1 (a), 0.2 (b), 0.3 (c) and 0.5 (d) Finally, to analyse the capabilities of the approach just described, for treating more realistic models, a MBS -FENE model composed of many dumbbells in a 2D physical space, is considered. This test is carried out in a simple shear flow with the following simulation conditions: b = 9, 1 du dy , 0.01 t , N = 5, M=2, Cn = 52, = 0.8 and being the initial condition given by Equation [14]. In all the simulations the flow kinematics is assumed known and uncoupled with the polymer presence.

The solutions at times t1 = 3 and t2 = 10 (at t2 the steady state is reached) are depicted in Figure 8. To simplify the solution post-processing, a smoothed distribution involving the connectors belonging to the same level (having the same index in the chain) is defined according to: From the computed results, we notice that symmetric levels (in our case levels 1 and 5, or 2 and 4) exhibit the same behaviour. Connectors are more stretched in the centre of the chain than on its ends. Figures 8(a) and 8(b) depicts 1 and 3 at time t1 and Figures 8(c) and 8(d) the same distributions at time t2. We can also notice the perfect symmetry of the computed distributions since the initial distributions respect symmetry. Figure 9 depicts the chain configurations at time t2 of some of the chains involved in the simulation. In this figure the origin of the first spring connector q 1 of all chains is located at the centre of the figure.

Conclusions

Despite of its simplicity, the deterministic particle technique just described seems to be robust and accurate when the chain model consists of a reduced number of connectors. However a number of difficulties persist, being the most important:

-the number of particles increases significantly with the dimension of the conformation space as was also found in the context stochastic simulations;

-the smoothing parameter has a strong impact on the computed results, even if it can be identified from the equilibrium conditions, when the distribution becomes highly localized an adaptation of this value is required to guaranty the solution accuracy. From some numerical experiments we have conclude that each particle must have a minimum number of neighbours that seems to be around 3;

-it is well known that in highly multidimensional spaces the use of a Euclidean distance results delicate, because in a high dimensional space the measure of the surface of a hyper-sphere becomes higher than the measure of its volume.

  [3] can be rewritten excluding Brownian effects which are retained in Equation[5]. Then, Equations [3] and [5] can be expressed, in a dimensionless form (seeBird et al., 1987, for more details), as:

.

  At the initial time the particles are distributed for verifying the equilibrium initial condition given by:

Figure 3 .

 3 Figure 3. Zoom of the distribution functions in the region close to the maximum moldeule stretching

Figure 4 .

 4 Figure 4. Histogram related to the particle position at t=1

Figure 6 .Figure 7 .

 67 Figure 6. Particles distribution related to the solution depicted in figure 5(right)

Figure 8 .

 8 Figure 8. MBS conformation in a simple shear flow

Figure 9 .

 9 Figure 9. Example of chain representation at time t2

  Ammar A., Chinesta F., A Particle Strategy for Solving the Fokker-Planck Equation Governing the Fiber Orientation Distribution in Steady Recirculating Flows Involving Short Fiber Suspensions, Lectures Notes on Computational Science and Engineering, Springer, 43, 2005, p. 1-16. Bird R.B., Crutiss C.F., Armstrong R.C., Hassager O., Dynamic of polymeric liquid, vol. 2, Kinetic Theory, John Wiley & Sons, 1987.


Jill-Jênn Vie

Fabrice Popineau
email: popineau@lri.fr

Yolaine Bourda
email: bourda@lri.fr

Éric Bruillard
email: eric.bruillard@ens-cachan.fr

Adaptive Testing using a General Diagnostic Model

In online learning platforms such as MOOCs, computerized assessment needs to be optimized in order to prevent boredom and dropout of learners. Indeed, they should spend as little time as possible in tests and still receive valuable feedback. It is actually possible to reduce the number of questions for the same accuracy with computerized adaptive testing (CAT): asking the next question according to the past performance of the examinee. CAT algorithms are divided in two categories: summative CATs, that measure the level of examinees, and formative CATs, that provide feedback to the examinees at the end of the test by specifying which knowledge components need further work. In this paper, we formalize the problem of test-size reduction by predicting student performance, and propose a new hybrid CAT algorithm GenMA based on the general diagnostic model, that is both summative and formative. Using real datasets, we compare our model to popular CAT models and show that GenMA achieves better accuracy while using fewer questions than the existing models.

Introduction

Computerized assessments are becoming increasingly popular. Meanwhile, the Obama administration has urged schools to make exams less onerous and more purposeful [START_REF] Zernike | Obama administration calls for limits on testing in schools[END_REF]. To reduce over-testing, we need to optimize the time spent on tests, asking only informative questions about the learners' ability or knowledge. This is the idea behind Computerized Adaptive Testing (CAT): selecting the next question to ask according to the previous answers of the examinee. As an example, 238,536 such adaptive tests have been administered by the Graduate Management Admission Council in 2012-2013 [START_REF] Council | Profile of GMAT Candidates -Executive Summary[END_REF] and adaptive assessment is getting more and more necessary in the current age of massive online open courses (MOOC), in order to minimize dropout.

Primarily, CATs have been relying on item response theory, that provides a framework to measure effectively scores called latent traits in order to rank students on a scale. The idea is to calibrate the difficulty of questions using a history of people having already taken the test. In 2001, The No Child Left Behind Act has called for more formative assessments, providing feedback to learners and teachers at the end of the test. Such formative assessments may detect students with cognitive disabilities or simply build a profile that specifies which knowledge components seem to be mastered and which ones do not. A straightforward application would be a personal assistant that asks a few questions, then highlights the points that need further work, and possibly suggests useful material for remediation. In 2003, to address this need, new CATs have been developed relying on cognitive diagnosis models, the most popular being the DINA model [START_REF] Huebner | An overview of recent developments in cognitive diagnostic computer adaptive assessments[END_REF] based on a q-matrix: a matrix that maps items (aka questions) to knowledge components involved in their resolution. Other cognitive models, less known, tend to unify scoring and formative assessments, but to date, they have not be used for adaptive testing [START_REF] Yan | Computerized multistage testing[END_REF].

In this paper, we formalize the problem of test-size reduction by predicting student performance (TeSR-PSP), inspired by [START_REF] Lan | Sparse factor analysis for learning and content analytics[END_REF] and present a new algorithm for CAT called GenMA, based on the general diagnostic model [START_REF] Davier | A general diagnostic model applied to language testing data[END_REF] that encompasses both the recovery of the latent knowledge components (KC) and for each KC, a degree of proficiency represented by a difficulty parameter.

To compare our algorithm for adaptive testing to the other ones mentioned above, we present an experimental protocol and execute it on real data. We show that GenMA outperforms existing models.

Our paper is organized as follows. In Section 1, we present the related work in CAT models. In Section 2, we formalize the problem of test-size reduction by predicting student performance and our new model, GenMA. In Section 3, we present an experimental protocol devised to compare these models, the real dataset used for evaluation, and our results. Finally, we discuss further work.

Related Work

Computerized Adaptive Testing and the problem of test-size reduction

In a non-adaptive test, every examinee receives the same set of questions. In an adaptive test, the next item asked by the system is chosen according to a certain criterion (the selection item rule), until the termination criterion holds, for example until a threshold over the parameters of the chosen model is guaranteed. Therefore, asking questions in an adaptive way means that the next question can be chosen as a function of the previous responses of the examinee.

In the problem of test-size reduction [START_REF] Lan | Sparse factor analysis for learning and content analytics[END_REF], one wants to reduce the number of questions asked as much as possible. Given a student model, we thus need to carefully choose the next question in order to still recover the model parameters. Formally, we want to decrease as much as possible the distance between the estimated and true user parameters after each question.

But in real data analysis, the true user parameters are unknown. For their evaluation, [START_REF] Lan | Sparse factor analysis for learning and content analytics[END_REF] replace the true user parameters with the estimated parameters they obtain after all questions have been asked, even if those estimated parameters do not fit the data at all.

In what follows, we will assume n learners take a test of a total of m questions. We assume the student data is dichotomous, which means every student either fails on succeeds over an item.

Item Response Theory (IRT)

The most simple model in item response theory for adaptive testing is the Rasch model, also known as 1-parameter logistic model. It models the behavior of a learner i ∈ {1, . . . , n} with a single parameter θ i ∈ R called ability, and models the item j ∈ {1, . . . , m} with a single parameter d j ∈ R called difficulty. The tendency for a learner to solve an item only depends on the difference between the difficulty and the ability:

P r("learner i answers item j") = Φ(θ i -d j )
where Φ : x → 1/(1 + e -x ) is the logistic function. Being a unidimensional model, the Rasch model alone is not suitable for finegrained feedback: it can only provide the level of the examinee at the end of the test. Still, it is really popular because of its simplicity, its stability and its sound mathematical framework [START_REF] Desmarais | A review of recent advances in learner and skill modeling in intelligent learning environments[END_REF][START_REF] Bergner | Model-based collaborative filtering analysis of student response data: Machinelearning item response theory[END_REF]. Also, [START_REF] Verhelst | Profile analysis: A closer look at the PISA 2000 reading data[END_REF] has showed that if the items are splitted into categories, the Rasch model is enough to provide to the examinee a useful deviation profile, specifying which category subscores were lower or higher than expected.

It is natural to extend the Rasch model to multidimensional abilities. In Multidimensional Item Response Theory (MIRT) [START_REF] Reckase | Multidimensional item response theory[END_REF], both learners and items are modelled by vectors of a certain dimension d, and the tendency for a learner to solve an item depends only on the dot product of those vectors. Thus, if learner i ∈ {1, . . . , n} is modelled by vector θ i and item j ∈ {1, . . . , m} is modelled by vector d j : P r("learner i answers item j") = Φ(θ i • d j ). Thus, a learner has greater chance to solve items correlated with its ability. Nevertheless, those richer models involve many more parameters, and have proven to be much harder to converge [START_REF] Desmarais | A review of recent advances in learner and skill modeling in intelligent learning environments[END_REF]. [START_REF] Xu | A simulation study to compare CAT strategies for cognitive diagnosis[END_REF] have used adaptive testing strategies applied to cognitive diagnosis (CD) models, notably the DINA model. These cognitive models rely on a specification of the knowledge components (KCs) involved in the resolution of the items proposed in the test, in the form of a q-matrix, which simply maps items to KCs: q ik is 1 if item i involves the KC k, 0 otherwise. Several algorithms have been proposed and compared for CD-CATs, using for example Kullback-Leibler divergence [START_REF] Huebner | An overview of recent developments in cognitive diagnostic computer adaptive assessments[END_REF][START_REF] Cheng | When cognitive diagnosis meets computerized adaptive testing: CD-CAT[END_REF].

Cognitive Diagnosis

If there are K KCs involved in a test, the learner can be modelled by a vector of K bits called state, specifying which KCs are mastered and which ones do not. Knowing the state of a learner, we can infer his performance over the different questions of the test. Slip and guess parameters capture careless errors. Throughout the assessment, a probability distribution over the 2 K states is maintained, and updated after each answer in order to fit the learner's behavior. In the particular case of the DINA model, the KCs involved in the resolution of an item are required to solve it. If the learner masters all KCs required, it still has a probability to slip over the question; if it lacks a KC, it still has a probability to guess correctly the answer.

Our Contribution

TeSR-PSP: Test-Size Reduction by Predicting Student Performance

In this paper, we propose a new problem called test-size reduction by predicting student performance: if we can ask only k questions in an adaptive way, which ones should we pick so as to predict the examinee's performance over the remaining questions of the test? Usually, adaptive tests keep going until a suitable confidence interval over the learner parameters is obtained. In our case, we want to specify in advance the maximal number of questions that will be asked to every student, in order to prevent boredom from the learner.

GenMA: Using the General Diagnostic Model for Adaptive Testing

[6] has proposed a unified model that takes many existing IRT models and cognitive models as special cases: the general diagnostic model for partial credit data: P r("learner i answers item j")

= Φ β i + K k=1 θ ik q jk d jk
where K is the number of KCs involved in the test, β i is the main ability of learner i, θ ik its ability for KC k, q jk is the (j, k) entry of the q-matrix which is 1 if KC k is involved in the resolution of item j, 0 otherwise, d jk the difficulty of item j over KC k. Please note that this model is similar to the MIRT model specified above, but only parameters that correspond to a nonzero entry in the q-matrix are taken into account.

To the best of our knowledge, this model has not been used in adaptive testing [START_REF] Yan | Computerized multistage testing[END_REF]. This is what we present in this paper: GenMA relies on a general diagnostic model, thus requires the specification of a q-matrix by an expert. The parameters d jk for every item j and KC k are calibrated using the history of answers from a test and the Metropolis-Hastings Robbins-Monro algorithm [START_REF] Chalmers | mirt: A multidimensional item response theory package for the R environment[END_REF][START_REF] Cai | Metropolis-Hastings Robbins-Monro algorithm for confirmatory item factor analysis[END_REF]. For the selection item rule of GenMA, we choose to maximize the Fisher information at each step, details of the implementation can be found in [START_REF] Chalmers | mirt: A multidimensional item response theory package for the R environment[END_REF]. The problem TeSR-PSP becomes: after k questions asked to a certain learner i, how to estimate its main ability β i and ability for each KC θ ik that can explain its behavior throughout the test?

In real tests, items usually rely on only few KCs, hence there are fewer parameters to estimate than in a general MIRT model, which explains why the convergence is easy to obtain for GenMA. We can thus use the general diagnostic model to create an adaptive test that makes best of possible worlds: providing feedback under the form of degrees of proficiency over several KCs at the end of test, represented by the vector θ i = (θ i1 , . . . , θ iK ), and being easy to converge. GenMA is both summative and formative, thus a hybrid model. Such feedback can be aggregated at various levels (e.g., from student, to class, to school, to city, to country) in order to enable decision-making [START_REF] Verhelst | Profile analysis: A closer look at the PISA 2000 reading data[END_REF][START_REF] Shute | Advances in the science of assessment[END_REF].

Evaluation

In this section, we detail the experimental protocol used to compare the following models for TeSR-PSP: the Rasch model, the DINA model and GenMA. For the sake of equality, we decide to define the same selection item rule for all models: all of them pick the question that maximizes Fisher information, which means the question of estimated probability closest to 1/2.

Experimental protocol

Our experimental protocol is based on double cross-validation. For each experiment, we need:

-a train student set, in order to calibrate the parameters of our model (for example, the difficulty of questions in the case of the Rasch model); -a test student set, which will take our adaptive test; -a validation question set V Q , which is used for training, but kept out of the adaptive tests, used only to evaluate the prediction of performance of the students from the test set.

To evaluate the score of a model for our problem, we first train it using the train student set. Then, for each student from the test set, we let the model pick a question, we reveal the student answer for this question, the model updates its parameters accordingly and outputs a probability of correctly answering the questions from the validation set, that we can evaluate using negative loglikelihood, hereby denoted as error : error(pred, truth) = q∈VQ truth q log pred q + (1truth q ) log(1pred q ). Then the model picks the second question, and so on. Thus, after k questions we can compute a prediction error over the validation question set for every model and every test student.

Real dataset: fraction subtraction

Tatsuoka's fraction subtraction dataset contains the dichotomous responses of 536 middle school students over 20 fraction subtraction test items. The corresponding q-matrix maps the 20 items to the following 8 knowledge components (KCs):

-convert a whole number to a fraction, -separate a whole number from a fraction, -simplify before subtracting, -find a common denominator, -borrow from whole number part, -column borrow to subtract the second numerator from the first, -subtract numerators, -reduce answers to simplest form.

The cross-validation was performed using a random split into training and test sets: the split was 5-fold over the students and 4-fold over the questions: each student set had a size of 20% while each validation question set had a size of 25% so 5 questions, therefore there were 20 experiments in total, of which the mean error was computed.

Implementation details

Our Rasch model implementation comes from the ltm R package. We made our custom implementation of the DINA model but the slip and guess calibration is held by the CDM package. GenMA is built upon the mirt package [START_REF] Chalmers | mirt: A multidimensional item response theory package for the R environment[END_REF].

Results

For each number of questions asked from 1 to 15, we plotted the mean error of each model (Rasch model, DINA model and GenMA) over the test student set in Figure 1, and as an insight easier to comprehend, the mean number of questions incorrectly guessed in Figure 2.

Figure 1 shows that 8 questions over 15 are enough for the Rasch model to converge on the fraction subtraction dataset. Figure 2 shows that no matter how many questions are asked, the Rasch and DINA models can't predict correctly more than 4 questions in average over 5 in the validation question set, while GenMA can achieve this accuracy with only 4 questions, then keeps on improving its predictions. The DINA model takes a long time to converge because the first questions require a single KC, therefore they do not bring a lot of information about the user state. But still, the simplest, unidimensional Rasch model performs surprisingly well compared to GenMA which is over 8 dimensions, one per KC of the q-matrix. In this paper, we formulated the problem of test-size reduction by predicting student performance, and presented our new adaptive testing algorithm GenMA to tackle it, based on the general diagnostic model. As this model is richer than other models as Rasch or DINA, it could be prone to overfitting: having more parameters, it may have a good score over train data but poor generalization over the test data. But we showed it actually achieves a better accuracy at predicting student performance, using fewer questions than the existing models on a real dataset.

The idea of a hybrid model combining several KCs and weights for each of them is not new: MIRT models can be seen this way, but there are many parameters to estimate, leading to convergence issues. [START_REF] Lan | Sparse factor analysis for learning and content analytics[END_REF] presented sparse factor analysis (SPARFA), a model that combine q-matrices and weights but their KCs are specified automatically, not by experts, thus it is not possible to provide a feedback at the end of the test.

In order to overcome the complexity of O(2 K ) of the DINA model, some knowledge representations such as Attribute Hierarchy Model [START_REF] Leighton | The attribute hierarchy method for cognitive assessment: A variation on Tatsuoka's rule-space approach[END_REF][START_REF] Rupp | Putting ECD into practice: The interplay of theory and data in evidence models within a digital learning environment[END_REF] or Knowledge Space Theory [START_REF] Doignon | Knowledge spaces[END_REF][START_REF] Lynch | Real world usage of an adaptive testing algorithm to uncover latent knowledge[END_REF] have been devised, relying on dependencies over KCs in the form of a directed acyclic graph. We would like to compare these approaches to GenMA. Also, richer models using ontologies [START_REF] Mandin | Basing learner modelling on an ontology of knowledge and skills[END_REF][START_REF] Kickmeier-Rust | Competence-based knowledge space theory. Measuring and Visualizing Learning in the Information-Rich Classroom[END_REF] are a promising direction of research.

Fig. 1 .

 1 Fig. 1. Comparing adaptive testing models. Evolution of error (negative log-likelihood) over the validation question set, after a certain number of questions have been asked.

Fig. 2 .

 2 Fig. 2. Comparing adaptive testing models. Evolution of the number of questions in the validation question set incorrectly predicted, after a certain number of questions have been asked.

Acknowledgements

This work is supported by the Paris-Saclay Institut de la Société Numérique funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02.


Nicolas Borie

Justine Falque

Product-Coproduct Prographs and Triangulations of the Sphere

Keywords: Tamari lattice, Catalan numbers, Prographs, Triangulation, Sphere

In this paper, we explain how the classical Catalan families of objects involving paths, tableaux, triangulations, parentheses configurations and more generalize canonically to a three-dimensional version. In particular, we present productcoproduct prographs as central objects explaining the combinatorics of the triangulations of the sphere. Then we expose a natural way to extend the Tamari lattice to the product-coproduct prographs.

Abstract. Dans ce papier, nous expliquons comment les familles d'objets combinatoires comptées par les nombres de Catalan, comme les chemins, tableaux, triangulations, parenthésages, ou autre, se généralisent dans le monde des objets comptés par les nombres de Catalan tridimensionnels. En particulier, nous présentons les prographes produit-coproduit comme centraux pour expliquer la combinatoire des triangulations de la sphère. Nous exposons alors une manière naturelle d'étendre le treillis de Tamari aux prographes produit-coproduit.

Introduction

This paper is about a generalization of the famous combinatorial problem defining Catalan numbers and involving parentheses configurations, binary trees and triangulations of the n-gon. This paper is a sequel of [START_REF] Borie | Three-dimensional Catalan numbers and product-coproduct prographs[END_REF] where product-coproduct prographs (PC prographs) were introduced to study the combinatorics of the three-dimensional Catalan world. The work presented in this paper follows the results obtained in FPSAC 2017 [START_REF] Borie | Three-dimensional Catalan numbers and product-coproduct prographs[END_REF], whose statements can be summarized in Table 1.

In this paper, we complete these first statements to the dual side of planar objects (triangulations of the sphere) and the combinatorics of the planar objects.

Mireille Bousquet-Mélou did expose in [START_REF] Bousquet-Mélou | Counting planar maps, coloured or uncoloured[END_REF] that the three-dimensional Catalan numbers count the number of triangulations of the bipolar sphere. By embedding PC prographs on the sphere and linking the only global output to the global input, we present PC prographs as Voronoï diagrams of the triangulations of the bipolar sphere. This implies that the original bijection from rectangular Young Tableaux to PC prographs can be extended to triangulations.

Integer sequences and formal power series Catalan numbers

Three-dimensional Catalan numbers (2n)! (n + 1)! • n! 2 • (3n)! n! • (n + 1)! • (n + 2)! First combinatorial classes of objects Two-rows standard Young tableaux Three-rows standard Young tableaux

• • • • • • 2n 1 • • • • • • • • • • • • 3n • • • • • • 1 • • • • • •

Realization as Operators rules in Algebras Parentheses configurations of an n-product

Ways to assemble-disassemble n-times each in an associative algebra in an associative and coassociative bialgebra

× • ∆ • × • ∆ (• × •) × • × • (× ⊗ Id) • (Id ⊗ ∆) • ∆ × • (Id ⊗ ×) • (Id ⊗ ∆) • ∆ • × (• × •) × • (× ⊗ Id) • (∆ ⊗ Id) • ∆ × • (Id ⊗ ×) • (∆ ⊗ Id) • ∆ Planar combinatorial objects Binary trees PC prographs
Table 1: Classical Catalan world and its three-dimensional counterpart.

By duality, since rotations in prographs correspond to triangle flipping, we define four rotations rules on PC-prographs. We explain that these rotations correspond exactly to all admissible triangulation flips on the triangulations side. Among these four rotations rules, two are extentions of the classical rotations balancing binary search trees and the other two are suggested by the geometry of the triangulations of the sphere. Finally, we present some structural properties of the set of PC prographs of size n extended with the four rotations viewed as oriented rewriting rules. We argue that, in some cases, PC prographs can be viewed as a gluing of two binary trees; endowing this subfamily with two of the four rotations makes it a lattice that is a product of the Tamari lattice by itself. However, although adding the other two rotations does generate all PC-prographs with faithful action of the rotations, this does not provide a lattice.

This paper is organised as follows. In the next Section, we first recall what a PC prograph is and describe the duality between PC prographs and rooted triangulations of the sphere. We exploit the duality to introduce and orient flips and rotations. In Section 4, we investigate structures that can be defined over the set of PC prographs.

Duality and triangulations

PC prographs are connected oriented planar graphs without cycles whose vertices are products or coproducts. Coproducts have one input and two outputs and products have two inputs and one output. Just one coproduct (resp. product) has its input (resp. output) unconnected, which is the input (resp. output) of the prograph (see Figure 2 for examples). Since each vertex has degree three, the dual graph consists of triangles. Now, embed a prograph on the sphere and connect its output to its input going through the exterior non visible side of the sphere. The prograph becomes a Voronoï diagram of some triangulation of the sphere. Since [START_REF] Bousquet-Mélou | Counting planar maps, coloured or uncoloured[END_REF]7] prove that bipolar oriented triangulations of the sphere are counted by the three-dimensional Catalan numbers, this procedure provides an alternative and explicit bijection.

For example, consider the following three-row and rectangular standard Young tableau: 7 8 11 12 3 4 9 10 1 2 5 6 . This tableau is stable by the Schützenberger involution (see [START_REF] Borie | Three-dimensional Catalan numbers and product-coproduct prographs[END_REF]6]) (rotate 180 • the tableau then complement its values sending i to 3n + 1i). Using the central bijection of [START_REF] Borie | Three-dimensional Catalan numbers and product-coproduct prographs[END_REF] from tableaux to prographs, we obtain the black prograph of Figure 1, with coproducts represented as circles and products as squares. Considering this prograph as a Voronoï diagram, we build the red oriented triangulation of the sphere with a single edge passing through its dark side. This red dashed edge (its midpoint could have been sent to infinity) splits the dark side of the sphere into two triangles: one containing the North pole of the sphere and the other one containing its South pole.

Remark 1. Switching the context from prographs to bipolar oriented triangulations of the sphere, the Schützenberger involution consists in rotating the sphere by 180 degrees, swapping the poles, and flipping the orientation of the edges of the triangulation. It is also equivalent to the antipodal map of the sphere with renamed poles since they have been swapped.

The reader can check that the triangulation obtained in Figure 1 is stable by the Schützenberger involution.

Theorem 1.

With PC prographs as an intermediate object, graph duality provides an explicit and constructive bijection from 3-row rectangular standard Young tableaux onto bipolar-oriented triangulations of the sphere. This bijection commutes with the Schützenberger involution.

Since there exists numerous algorithms for exhaustive enumeration and (possibly random) generation of Young tableaux, this bijection using PC prographs enables a lot of calculations and computations on the triangulations side. For example, we believe that the easiest way to check if two oriented triangulations of the bipolar sphere are homotopic is to build their associated prographs with graph duality and check whether they are associated with the same Young tableau.

Remark 2. The world of prographs meets here the world of oriented planar maps. Indeed, prographs are planar assemblies of operators with identified inputs and outputs; if one puts reasonable operators on the vertices of a planar oriented map, with respect to the number of ongoing and outgoing edges for each vertex, one builds a prograph. Since, for PC prographs, a single entry means the vertex is a coproduct and two entries mean it is a product, there is no choice to be made and everyone works on the same objects.

Flips and rotations

The Tamari lattice has numerous different realizations : using parentheses configurations, binary trees, ordered forests, triangulations and more. Our first attempts to uncover a structure on PC prographs were by focusing on the prographs and Young tableaux side. Unfortunately, there were too many choices when trying to establish a set of transforms and give them orientations, thereby exhibiting a partial order (or more).

On the other hand, the geometric approach revealed the needed rigidity. When one considers triangulations, the reasonable flip consists in focusing on an edge, check that it is the border of two triangles; if the union of these two triangles forms a quadrilateral and the edge is a diagonal of it, then one flips this edge by removing it and drawing the other diagonal instead. Using this rule alone and by exhaustion of cases, we now show that there is only one reasonable way to consider flips and rotations.

By graph duality, each triangle contains a single operator, a product or a coproduct. Pairs of neighboring triangles (sharing at least an edge) can be enumerated by the entries inside the associated Young tableau since, if triangles are neighbors, there exists an edge between both involed operators. We forget 1 which virtually connects the northernmost product with the southernmost coproduct via the dark side of the sphere. This pointed edge in PC prographs correspond to a pointed edge in triangulations that we will not flip (it plays the same role as the first side of the n-gon in the classical Catalan world and Tamari lattice.) In PC prographs, there are three types of outputs : outputs of products, left outputs of coproducts and right outputs of coproducts. Similarly, there are three types of inputs : inputs of coproducts, left inputs of products and right inputs of products. Since edges in PC prographs connect an output of operator to an input, we should consider nine cases.

The only choice we made is the orientation of operations. Here, types I, V and IX are reduced by choice, and, since we want to be coherent with the Schützenberger involution, the only other possible choice is to reverse all orientations simultaneously. Reversing all rules, types III, IV and VI would have been the reduced edges. This corresponds to exchanging the left and right everywhere and flipping the east and west on the sphere.

The reader can check exhaustively that a flip of types II and VI would create a loop inside the prograph, whatever the choice of orientation. Therefore flipping these two types of edges is technically impossible. It so happens that these edges correspond exactly to pathological organizations of triangles on the sphere. For example, on the sphere, there exist valid triangulations in which two neighboring triangles have two edges in common. The union of two such triangles does not form a quadrilateral and no flip is possible.

Theorem 2. The four rotation rules in PC prographs shown in Figure 5 correspond to all valid edge flippings in rooted oriented triangulations of the sphere.

The reader can check that the first two rotations are directly derived from the classical rotation rule in binary trees.

The four rotation rules, without their orientations, are completely rigid as far as geometry is concerned. Anyway, the reader could argue that choices of orientation were made. Since we still want to be compatible with the Schützenberger involution, the choice of minimum and maximum must be reversed for product and coproduct. But exactly as the classical Tamari lattice implies a choice (one can set either the left comb tree or the right comb tree to be the maximum), there is an equivalent choice with PC prographs. This will result in the same kind of symmetries for the classical Tamari lattice and our new structure on PC prographs.

Poset structure over PC prographs

In this section, we investigate what type of structure can be placed over the set of PC prographs.

Gluing two trees at their canopies

The orientation in rooted triangulations of the sphere come directly from orientations of the classical case, the triangulations of the rooted n-gon. Choosing a root side in the n-gon amounts to set an orientation to a chosen side (ongoing or outgoing, the flow being perpendicular to this side) and set the reverse orientation to all other sides. So there are only two isomorphic choices. On the one hand, one can set the root edge to be outgoing. The other sides are thus ongoing, hence inputs, and the dual of the triangulation ends up being a rooted tree of products. On the other hand, one can set the root edge to be ongoing, all other sides become outputs and the dual becomes a rooted tree of coproducts. Now, consider two trees with the same number of nodes. Build their associated triangulations of the n-gon and set opposite orientations on each of them (as illustrated in Figure 7). By deformation, we can imagine these two n-gons are two hemispheres that can be glued along their equator, the root sides being glued together. This way, considering the underlying trees, the n -1 outputs of the tree of coproducts end up connected to the n -1 inputs of the tree of products. Call the product of the outgoing root side the North and the coproduct of the ongoing root side the South, and set them at the top and bottom by convention, respectively. This operation produces a PC prograph having its products connected up to the North point and its coproducts connected up from the South. North and South are also connected by the only upside-down edge of the equator, on the dark side of the sphere. We obtain a figure similar to Figure 1.

An immediate property is that PC prographs having coproducts all connected from the South and products all connected up to the North are exactly PC prographs avoiding a pattern : the output of a product grafted as input of a coproduct. Theorem 3. The subset of PC prographs avoiding edges of types VII (connecting a product output to a coproduct input) endowed with the two classical flip types forms a lattice which is exactly the product lattice of the classical Tamari lattice with itself.

PC prographs whose coproducts show up first

In Section 5 of [START_REF] Borie | Three-dimensional Catalan numbers and product-coproduct prographs[END_REF], we did expose a process labelling the operators : from bottom to top and from left to right for operators whose entries lean on already labelled operators. This natural process (called Boriefication in [4]) is recalled in Figure 8. Among all PC prographs, a small portion produce labels 1, 2, . . . , n to coproducts and labels n + 1, n + 2, . . . , 2n to products. In some sense, the depth-left-first transversal of operators enumerate all coproducts before the first product can go out. Regarding standard tableaux, it is equivalent to check that the last entry of the first row (indexing the last coproduct) is smaller than the first entry of the last row (indexing the first product to have its two entries labelled). This subfamily of PC prographs is counted by Sequence A274969 of OEIS [START_REF] Sloane | The on-line encyclopedia of integer sequences[END_REF]. It is also a subfamily of that described in Section 4.1. Such prographs can be obtained as a gluing of two binary trees with a reattachment condition.

Labelling nodes of binary trees T with the depth-left-first search algorithm, we can refine the enumeration of binary trees by counting leaves to the left of the last node. Denoting by size(T) the number of nodes in T and by dg(T) (disabled grafting sites) the number of leaves to the left of the last node, we define Cat i (q) to be the following formal sum over all binary trees: Cat i (q) = ∑ size(T)=i q dg(T) .

(4.1)

Cat i (q) deploys one very famous Catalan triangle. Here are the first values.

Cat 0 (q) = 1 Cat 1 (q) = 1 Cat 2 (q) = 1 + q Cat 3 (q) = 1 + 2q + 2q 2 Cat 4 (q) = 1 + 3q + 5q 2 + 5q 3 Cat 5 (q) = 1 + 4q + 9q 2 + 14q 3 + 14q 4 Cat 6 (q) = 1 + 5q + 14q 2 + 28q 3 + 42q 4 + 42q 5

It is known that the coefficients of this triangle satisfy

C(n, k) = n + k k - n + k k -1 . Proposition 1.
∑ d 1 +d 2 n+1 C(n, d 1 )C(n, d 2 ) = 3n n -2 3n n -1 + 3n n -2
This formula is more complicated than the compact combination of three binomial due to Janis Stipins A274969 of OEIS [START_REF] Sloane | The on-line encyclopedia of integer sequences[END_REF]. Anyway, it provides an alternative way to compute OEIS A274969 : 1, 1, 4, 21, 121, 728, 4488, ... Furthermore, it exposes an alternative description of pushall stack words of length 3n [START_REF] Pierrot | 2-stack pushall sortable permutations[END_REF] as pairs of binary trees. The first values are the sum of the first coefficients of the square of Cat i (q). Cat 0 (q) 2 mod q = 1 = |q=1 1 Cat 1 (q) 2 mod q 2 = 1 = |q=1 1 Cat 2 (q) 2 mod q 3 = 1 + 2q + q 2 = |q=1 4 Cat 3 (q) 2 mod q 4 = 1 + 4q + 8q 2 + 8q 3 = |q=1 21 Cat 4 (q) 2 mod q 5 = 1 + 6q + 19q 2 + 40q 3 + 55q 4 = |q=1 121 Cat 5 (q) 2 mod q 6 = 1 + 8q + 34q 2 + 100q 3 + 221q 4 + 364q 5 = |q=1 728 Cat 6 (q) 2 mod q 7 = 1 + 10q + 53q 2 + 196q 3 + 560q 4 + 1288q 5 + 2380q 6 = |q=1 4488 4.3 A poset that is not a lattice on PC prograph Proposition 2. The set of all PC prographs with the four rotation rules does not form a lattice.

The generalization of the classical Catalan world stops here. The smallest pair of uncomparable elements are prographs of size 3. There are 3 elements higher in the poset than both prographs but no join can be defined (there are two uncomparable minima among these 3 higher elements). [4] Christophe Cordero. "Explorations combinatoires des structures arborescentes et libres". Theses. Université Paris-Est, Dec. 2019. url: https://tel.archives-ouvertes.fr/tel-0301

[5] The Sage Developers. Sage Mathematics Software (Version 6.4). http://www.sagemath.org. 2014.

[6] Justine Falque. "A bijection between weighted Dyck paths and 1234-avoiding alternating permutations". In: Sém. Lothar. Combin. 85B (2021), Art. 71, 13.

[

Figure 1 :

 1 Figure 1: A self-dual prograph (black) and its bipolar-oriented triangulation (red).

Figure 2 :Figure 3 :

 23 Figure 2:The five prographs with two coproducts, two products (in black), and their associated triangulations (in red).

Figure 4 :

 4 Figure 4: The 9 types of edges and their reductions via oriented rotations

Figure 5 :

 5 Figure 5: The four local rotations over PC prographs shown with their dual.

Figure 6 :

 6 Figure 6: The poset (lattice for this size) of PC prographs of size 2.

Figure 7 :

 7 Figure 7: Build the sphere and a PC prograph from two triangulations of the n-gon.

Figure 8 :

 8 Figure 8: Labelling operators of a prograph and its image by the Schützenberger involution.

Figure 9

 9 Figure9displays the poset over the 42 PC prographs of size 3 (as standard tableaux). This research was driven by computer exploration using the open-source mathematical software Sage [5] and its algebraic combinatorics features developed by the Sage-Combinat community[START_REF]The Sage-Combinat community. Sage-Combinat: enhancing Sage as a toolbox for computer exploration in algebraic combinatorics[END_REF].

Figure 9 :

 9 Figure 9: The poset formed by the 42 PC prographs of size 3.

  PC prographs of size n whose operators are labelled by indices 1, 2, . . . , n for coproducts and indices n + 1, n + 2, . . . , 2n for products are counted by pairs of binary trees of size n with d 1 and d 2 disabled grafting sites and d 1 + d 2 n + 1. Equivalently, we have


Member, IEEE Abi Waqas
email: abi.waqas@faculty.muet.edu.pk

Senior Member, IEEE Daniele Paolo Manfredi
email: paolo.manfredi@polito.it

Member, IEEE. Melati

P Manfredi

D Melati

Performance Variability Analysis of Photonic Circuits with Many Correlated Parameters

Keywords: Correlated manufacturing variability, Karhunen-Loève transform, performance prediction, photonic devices, polynomial chaos, principal component analysis, process variations, silicon photonics, uncertainty quantification. I

We propose a method to analyze the performance variability caused by fabrication uncertainty in photonic circuits with a large number of correlated parameters. By combining a sparse polynomial chaos expansion model with dimensionality reduction in the form of Karhunen-Loève transform and principal component analysis, we demonstrate the stochastic analysis of the transfer function of cascaded Mach-Zehnder interferometers with up to 38 correlated uncertain parameters.

Alternatively, stochastic spectral methods based on polynomial chaos expansion (PCE) have emerged as a promising alternative, showing significant speedup over Monte Carlo in many engineering problems [START_REF] Eldred | Recent advances in non-intrusive polynomial chaos and stochastic collocation methods for uncertainty analysis and design[END_REF]. PCE has been also proposed for the variability analysis of photonic circuits [START_REF] Weng | Stochastic simulation and robust design optimization of integrated photonic filters[END_REF], [START_REF] Waqas | Stochastic process design kits for photonic circuits based on polynomial chaos augmented macro-modelling[END_REF], [START_REF] Xing | Stochastic collocation for device-level variability analysis in integrated photonics[END_REF], [START_REF] Waqas | Efficient variability analysis of photonic circuits by stochastic parametric building blocks[END_REF], [START_REF] Cao | Hierarchical polynomial chaos for variation analysis of silicon photonics microresonators[END_REF], [START_REF] Cui | Chance-constrained and yield-aware optimization of photonic ICs with non-gaussian correlated process variations[END_REF]. While classical PCE assumes random parameters to be mutually independent, photonic devices exhibit spatial correlation over rather long distances of few millimeters that must hence be taken into account to properly analyse the impact of fabrication uncertainty [START_REF] Pond | Predicting the yield of photonic integrated circuits using statistical compact modeling[END_REF], [START_REF] Yang | Phase coherence length in silicon photonic platform[END_REF]. Some PCE implementations proposed for photonics problems take into account parameters correlation [START_REF] Xing | Stochastic collocation for device-level variability analysis in integrated photonics[END_REF], [START_REF] Cui | Chance-constrained and yield-aware optimization of photonic ICs with non-gaussian correlated process variations[END_REF], but they either consider only a small set of random parameters or require the complex calculation of specialized polynomial basis. Indeed, one of the major drawbacks of PCE is that its computational cost grows exponentially with the number of stochastic variables, thus eventually becoming less efficient than Monte Carlo.

In this work, we tackle these issues by proposing an effective PCE implementation suitable for photonic circuits that can handle a large number of correlated uncertain parameters. The proposed method leverages the combination of several elements. Specifically, Karhunen-Loève transform (KLT) [START_REF] Loève | Probability Theory[END_REF] is applied to decorrelate the physical parameters and reduce their number, thereby also reducing the complexity of the model to be trained. Principal component analysis (PCA) is then used to compress the training data and reduce the output quantities to be modeled. Therefore, the entire frequencydependent transfer function is reduced to a small number of "principal components", for which an individual model is trained. An adaptive sparse PCE model [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF] is built for the principal components, from which the complete transfer function is recovered via the inverse PCA transformation. The sparse implementation is particularly suited for training highdimensional and/or high-order PCEs with a small amount of training data. Two test cases, i.e., a fifth-order and a seventhorder photonic filter realized with cascaded Mach-Zehnder interferometers, are used to demonstrate the proposed method. Filters have up to 38 correlated uncertain parameters and the variability analysis of their entire frequency-domain transfer function is performed more than 30 times faster than classical Monte Carlo analysis.

II. PROPOSED METHOD

We consider a generic system that depends on a set of d uncertain input parameters x = (x 1 , . . . , x d ) T ∈ R d . We represent the system as

y = M(x), (1)
where M : R d → R is an implicit function that maps a given configuration of the input parameters x to the corresponding output y of the system. For the moment, we assume that the output of ( 1) is a scalar. We will later relax this assumption. In the context of this paper, x could be for example a collection of critical physical parameters of an optical device, e.g., the waveguide width or the gap of a directional coupler, and y an output of interest such as the transfer function at a given wavelength, the splitting ratio of a power splitter, or the bandwidth of a filter. The model can be equally applied to isolated components or entire photonic circuits.

A. Classical PCE

At this stage, we further assume that the parameters x are statistically independent. The classical PCE framework seeks an approximation of (1) in the form of

y ≈ M(x) = α∈A c α ϕ α (x) (2)
where A is a set of multi-indices defined as [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF]]

A = A d,p q = α = (α 1 , . . . , α d ) ∈ N d : α q ≤ p , (3)
ϕ α are multivariate orthonormal polynomials, and c α are the pertinent coefficients. The polynomial basis functions are constructed as the product

ϕ α (x) = d j=1 φ αj (x j ), (4)
where φ αj is a univariate polynomial in the variable x j , with deg(φ αj ) = α j . Therefore, the multi-indices define the order of the basis function in each variable. For example, the polynomial ϕ (2,0,1) is quadratic in x 1 , constant in x 2 , and linear in x 3 .

The univariate polynomials φ αj in (4) are orthonormal based on the inner product

φ αj , φ αm = R φ αj (x j )φ αm (x j )w(x j )dx j = δ αj ,αm , (5)
where w(x j ) is the probability density function of x j . This ensures also the orthogonality of the multivariate polynomials ϕ α , as well as an exponential decay of the error norm

M -M 2 2 = R d M(x) -M(x) 2 w(x)dx (6)
for functions with finite variance [START_REF] Xiu | The wiener-askey polynomial chaos for stochastic differential equations[END_REF], where w(x) is the joint probability density function of the variables x, as the PCE order p is increased. It should be noted that the error ( 6) is defined in statistical terms: a larger deviation is tolerated for values of x with lower probability. A useful property of PCEs is that the first two statistical moments (i.e., the mean and the variance) are available in closed form as a function of the coefficients, i.e.,

E{y} ≈ E{ M} = R d M(x)w(x)dx = c 0 , (7)
with 0 denoting the null element of N d , and

Var{y} ≈ Var{ M} = R d M(x) -c 0 2 w(x)dx = α∈A\0 |c α | 2 (8)
Several choices are commonly used for choosing the indices in (3) based on the parameter q. The most popular is certainly q = 1, leading to the so-called "total-degree" truncation, with a number of terms in (2) (i.e., the cardinality of A) given by

|A| = |A d,p 1 | = (p + d)! p!d! . (9)
Another common choice is q = ∞, which however leads to excessively large expansions. Instead, setting q < 1 provides increasingly sparser expansions that neglect higherorder interactions between input parameters, which is often acceptable [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF].

B. Calculation of PCE Coefficients

There exist several strategies for the calculation of the PCE coefficients [START_REF] Xiu | Fast numerical methods for stochastic computations: A review[END_REF]. The most popular are the stochastic Galerkin method, which is intrusive and requires to modify the system of governing equations, numerical integration techniques that approximate projection integrals, and least-square regression or interpolation approaches. Here, we focus on the third ones, which are the most straightforward, versatile, and scalable.

Consider a set {x } L =1 of "training samples" of the input parameters, and a vector of corresponding observations y = (y 1 , . . . , y L ) T ∈ R L , with y = M(x ), ∀ = 1, . . . , L. The least-square regression minimizes the norm of the residual between the model prediction at the training samples and the observations. The solution is cast as

c * = arg min c Ψc -y = Ψ + y, (10)
where c * ∈ R K is a vector collecting all PCE coefficients, Ψ ∈ R L×K is a matrix with entries Ψ ,α = ϕ α (x ), and Ψ + = (Ψ T Ψ) -1 Ψ T is its Moore-Penrose pseudo-inverse, which generalizes the matrix inverse for non-square matrices.

For an acceptable solution, the regression problem needs to be overdetermined, i.e., L > K (typically, at least two to three times larger). Plain random sampling, Latin hypercube sampling (LHS), or quasi-MC low-discrepancy sequences [START_REF] Singhee | Why quasi-monte carlo is better than monte carlo or latin hypercube sampling for statistical circuit analysis[END_REF] can be used to generate the training samples, with an impact on the convergence performance.

C. Sparse PCE

The standard approach rapidly becomes computationally intractable when the number of input parameters increases. For example, with an expansion order as low as p = 2, d = 30 variables, and q = 1, the number of basis functions is K = 496. Therefore, a number of training samples in the order of 10 3 is needed, with a computational effort that becomes comparable to a MC analysis. If order p = 4 is required, the number of training samples becomes of the order of 10 5 , while remaining of the order of a few thousands even if a sparser expansion with q = 0.7 is adopted.

Regardless of the initial sparsity of set A d,p q , typically most of the PCE coefficients further turns out to be negligible as a consequence of the "sparsity-of-effects principle" [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF]. Therefore, sparse approaches were developed to identify only nonnegligible coefficients, with a reduced requirement in terms of training samples. This allows extending the PCE framework to problems characterized by a large number of uncertain input parameters and/or requiring high-order expansions.

In this paper, we make use of the UQLab toolbox [START_REF] Marelli | Uqlab: A framework for uncertainty quantification in matlab[END_REF] to calculate sparse PCEs. The toolbox implements the adaptive method in [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF], which is based on least-angle regression. Given the dimensionality d, the truncation parameter q, and a vector of possible p-values, the algorithm identifies the optimal expansion order p as well as the coefficients of a sparser subset of basis functions à ⊂ A d,p q , with | Ã| |A d,p q |. We refer the Reader to the toolbox documentation for further information.

D. Correlated Variables: Karhunen-Loève Transform

The method so far assumes that the uncertain input parameters are mutually independent. In the most general case of parameters with arbitrary correlation, customized basis functions can be numerically computed [START_REF] Cui | Stochastic collocation with non-gaussian correlated process variations: Theory, algorithms and applications[END_REF]. However, for the special case of Gaussian correlated variables, KLT can be effectively used to decorrelate the variables [START_REF] Loève | Probability Theory[END_REF].

Assume that x are correlated Gaussian variables with mean vector µ ∈ R d and covariance matrix Σ ∈ R d×d . The variables x can be expressed in terms of a set of independent and standard (i.e., with zero mean and unit variance) Gaussian

variables ξ = (ξ 1 , . . . , ξ d ) T ∈ R d as x = µ + d k=1 λ k ξ k p k (11)
where λ k and p k are the eigenvalues and corresponding eigenvectors of the covariance matrix Σ. It should be noted that a trivial application of ( 11) is normally used to renormalize independent uncertain parameters to a common support, so that a standard set of univariate basis functions can be used for each variable. With the transformation 11, the PCE framework introduced so far can be seamlessly applied in terms of the variables ξ. Furthermore, for highly-correlated variables, the eigenvalues decay very fast. Hence, the KLT (11) can be truncated in such a way that the d variables x are expressed in terms of a reduced set of d independent random variables [START_REF] Ochoa | Random-space dimensionality reduction for expedient yield estimation of passive microwave circuits[END_REF], i.e.,

x ≈ µ + d k=1 λ k ξ k p k ( 12
)
where the eigenvalues are now assumed to be sorted in decreasing order and d is such that

λ k λ 1 ≤ ε KLT , (13)
for k > d and a prescribed threshold ε KLT . This dimensionality reduction allows alleviating the computational complexity of the PCE model.

E. Multiple Outputs: Principal Component Analysis

Another difficulty arises when the output of system (1) is not a scalar. This is the case, for instance, when a transfer functions is computed at different wavelengths within a given frequency band. The solution of the standard regression approach ( 10) is readily applied in a vectorized manner to the output vector thanks to the fact that the basis functions, and hence the regression matrix Ψ, are the same for each component. On the contrary, sparse methods, such as that considered in this work, identify the basis functions based on the output data, and therefore the subset à potentially differs for each component of the vector of the outputs. One trivial solution is to apply the method individually to each output, which however is computationally inefficient when their number is large.

An effective strategy is to compress, using PCA, the output variables into a reduced set of components, for which the training of individual PCE models becomes feasible [START_REF] Manfredi | A data compression strategy for the efficient uncertainty quantification of time-domain circuit responses[END_REF]. Let us assume that a set of L training observations are available for M outputs of interest, and stored into a matrix Y ∈ R M ×L . Using an "economy-size" singular value decomposition (SVD), the matrix dataset is expressed as [START_REF] Yaghoubi | Sparse polynomial chaos expansions of frequency response functions using stochastic frequency transformation[END_REF]

Y = ȳ + U SV T = ȳ + U Z ( 14
)
where ȳ ∈ R M is the columnwise mean of Y over the training observations, U ∈ R M ×L , V ∈ R L×L , and S ∈ R L×L is a diagonal matrix collecting the singular values {σ i } L i=1 in decreasing order.

Since the singular values typically decay very fast due to the interdependency between the outputs, the SVD is truncated to retain only ñ "principal components", leading to

Y ≈ ȳ + ñ i=1 u i z i ( 15
)
where {u i } M i=1 ∈ R M are the first ñ left-singular vectors of Yȳ, whereas {z i } ñ i=1 ∈ R L are in fact observations of the principal components. Similarly to the truncated KLT [START_REF] Cao | Hierarchical polynomial chaos for variation analysis of silicon photonics microresonators[END_REF],

ñ is selected such that σ i σ 1 ≤ ε PCA (16)
for i > ñ. The KLT and the PCA are in fact equivalent, since the latter can be cast as the eigenvalue decomposition of the experimental covariance matrix of the dataset Y [START_REF] Yaghoubi | Sparse polynomial chaos expansions of frequency response functions using stochastic frequency transformation[END_REF]. However, the purpose of KLT is to decorrelate and reduce the number of input parameters, whereas PCA is used to compress output data. In all simulations, we set ε KLT = ε PCA = 10 -3 . The observations {z i } ñ i=1 are used to train individual PCE models for the principal components, which is feasible because typically ñ M . Once a model is available for the principal components, it is used to generate predictions thereof, while predictions of the original outputs are recovered via [START_REF] Loève | Probability Theory[END_REF].

III. APPLICATION AND NUMERICAL RESULTS

The method proposed in Sec. II can be applied to any photonic circuit design, including for example resonant structures, since no assumption is done in its formulation on the type of system response [START_REF] Smit | An introduction to InP-based generic integration technology[END_REF]. Here, a fifth-order and a seventhorder photonic filter realized with cascaded Mach-Zehnder interferometers are chosen as examples to demonstrate the performance of the proposed approach for sets of random variables of different sizes. In the following, we consider as a reference the typical parameter values for the silicon-oninsulator technology. The complex frequency response (transfer function) of the photonic filters is the stochastic process under investigation. The fabrication uncertainty is represented as correlated and Gaussian-distributed random variables.

A. Fifth-Order Unbalanced Mach-Zehnder Filter

The fifth-order Mach-Zehnder filter is assembled using two standard elementary building blocks, i.e., waveguides and directional couplers, as schematically shown in Fig. 1. The nominal design of the filter is obtained with the synthesis technique described in [START_REF] Madsen | Optical Filter Design and Analysis: A Signal Processing Approach Optical Filter Design and Analysis: A Signal Processing Approach[END_REF]. The nominal 3-dB bandwidth is BW 0 = 39 GHz and the coupling coefficients for the six directional couplers are K 1 = K 6 = 0.993, K 2 = K 5 = 0.03, and K 3 = K 4 = 0.421, with an in-band isolation larger than 30 dB. The dependence of the coupling constant of the directional couplers on waveguide gap, width, and thickness is first determined using coupled mode theory and preliminary electromagnetic simulations. For the nominal case, a waveguide width of 500 nm, thickness of 220 nm, and gap of 300 nm allow obtaining the required coupling coefficients K 1 -K 6 with the coupling lengths Lc 1 = Lc 6 = 26.38 µm, Lc 2 = Lc 5 = 3.09 µm, and Lc 3 = Lc 4 = 12.53 µm. The effect of the transition regions in the couplers is neglected. For the waveguide building block, nominal width and thickness are again 500 nm and 220 nm, respectively, corresponding to an effective index and group index of 2.44 and 4.18. All the five stages of the filter have the same unbalance length of 681 µm, corresponding to a free spectral range of 100 GHz. Figure 2 shows, with bold lines, the ideal transfer function of the Mach-Zehnder filter at the bar and the cross ports, respectively.

Next, different sources of fabrication uncertainties are considered. For the waveguides in the two branches of each filter stage we assume an uncertainty only in the waveguide width, resulting in variability of the effective and group indices of the waveguide mode. The main effect of this variability will be a wavelength shift of the response of each stage. For the directional couplers, we consider a larger number of uncertainties, i.e., on the width and thickness of the waveguides and on the separating gap. Their collective effect is expected to mainly influence the resulting power coupling coefficients. [START_REF] Waqas | Sensitivity analysis and uncertainty mitigation of photonic integrated circuits[END_REF]. The vector x of the uncertain input parameters is then split into two components, x 1 and x 2 . The former collects parameters W gi (waveguide width of the i-th coupler), g i (coupling gap of the i-th coupler), W Uj (width of the upper arm of the j-th Organizing the parameters in two different groups allows defining a rather complex and realistic correlation scenario. The parameters in either of the two vectors x 1 and x 2 are assumed to be mutually correlated among themselves, while no correlation is introduced between parameters belonging to different sets. This reflects the fact that the source of uncertainty is different for the two sets of parameters, mainly originating from lithography and etching for the first group, and from wafer fabrication for the second.

The pertinent correlation coefficients are calculated by considering the correlation function [START_REF] Wagner | Monte carlo simulation of electromagnetic scattering from two-dimensional random rough surfaces[END_REF]

g(z) = e -z 2 /2l 2
, where l is the correlation length and z is the distance between the circuit components along the horizontal direction. As an example, the vertical dashed lines in Fig. 1 highlight the distance between the first and the second directional coupler, which is used for the definition of z in the calculation of the correlation between the parameters of these two components. A realistic correlation length of l = 4.5 mm is considered for both parameter groups [START_REF] Yang | Phase coherence length in silicon photonic platform[END_REF]. Since the distance between the upper and lower waveguide arms is relatively small in an unbalanced Mach-Zehnder filter, a correlation coefficient of one is considered between their parameters. This means that, in the simulation, the same variation is assigned to the upper and lower waveguides.

B. Seventh-Order Unbalanced Mach-Zehnder Filter

The seventh-order Mach-Zehnder filter is assembled in the same way as described above, with seven delay sections and eight directional couplers. The nominal design exhibits a nominal 3-dB bandwidth of BW 0 = 29 GHz. The coupling coefficients of the eight directional couplers are K 1 = K 8 = 0.9914, K 2 = K 7 = 0.04, K 3 = K 6 = 0.0690, and K 4 = K 5 = 0.3489, with an in-band isolation larger than 20 dB. Width and thickness of the waveguides are the same as for the fifth-order filter. The nominal gap of the directional couplers is again 300 nm, leading to coupling lengths Lc 1 = Lc 8 = 26.22 µm, Lc 2 = Lc 7 = 0.35 µm, Lc 3 = Lc 6 = 4.71 µm, and Lc 4 = Lc 5 = 11.21 µm. All the seven stages have again the same unbalance length of 681 µm, corresponding to a free spectral range of 100 GHz. For each building block, the same geometrical parameters as for the previous test case are considered to be uncertain, with the same probability distribution and correlation. Owing to the increased number of stages, the total number of uncertain parameters is now d = 38.

C. Simulation Results

Circuit analysis is performed using circuit-level simulations. In each filter stage, waveguide sections (including both upper and lower waveguide arms) and directional couplers are building blocks with 2 input and 2 output ports that are represented by means of 2 × 2 transfer matrices with complex coefficients [START_REF] Madsen | Optical Filter Design and Analysis: A Signal Processing Approach Optical Filter Design and Analysis: A Signal Processing Approach[END_REF]. The covariance matrix of the input parameters is first computed as described in Section III-A, taking into account the relative distance of each building block within the circuit. For each simulation run, the value of the parameters is then determined according to their mean and the computed covariance matrix, and used to calculate the coefficients of the transfer matrix of each building block. Lastly, transfer matrices are multiplied according to the building block connections to obtain the complex transfer function of the circuit in the spectral domain. In each simulation, 1000 wavelength samples are used.

For both test circuits, a Monte Carlo analysis is performed running 10000 simulations and it is used as a reference. As an example, Fig. 3 (a) and (c) report, with thin gray lines, the absolute value of some of the computed transfer functions (transmission) of the fifth-order filter at the bar and cross ports, respectively. The black dashed line marks λ = 1521.8 nm, which is the designed central wavelength of the filter. As can be seen, the transmission is severely affected, especially with large fluctuations of the 3-dB pass-band as well as a reduced isolation both in band (cross port) and out of band (bar port). Remarkably, these large fluctuations results from relatively small variations of the geometrical parameters. Indeeed, a variation of 10 pm in the waveguide width corresponds to an effective index variation of 8 • 10 -5 , whereas a coupling gap variation of 1 nm corresponds to a power coupling coefficient variation of 5 • 10 -3 . In the two figures, the mean (blue solid line, computed at each wavelength as the absolute value of the mean of the complex transfer function) and standard deviation (blue dashed line, again computed using the complex transfer function) of the Monte Carlo samples are reported as well. To further investigate the effect of variability, Fig. 3(b) and(d) show the probability density function of the transmission at the bar and cross port, respectively, obtained from the Monte Carlo samples at λ = 1521.8 nm (blue bars). While the transmission at the bar port is confirmed to be almost unaffected by the variability, the isolation at the cross port varies across a range of more than 30 dB, with -20 dB as the most probable value compared to the -32 dB of the ideal design without variability. In order to test the PCE-based methods, we first apply KLT to the 28 input correlated variables and we train a classical PCE model to compare its performance with our proposed method. By considering a threshold of ε KLT = 10 -3 for the truncation, the original correlated variables are reduced to a set of d = 5 independent variables. A different PCE model needs to be trained for each of the 1000 wavelength samples of the complex transfer function of the filter at the bar and cross port. For each model, we select an order p = 6 and q = 1, resulting in K = 462 coefficients, and we use L = 3K = 1386 samples for the training, drawn according to a LHS design. It should be noted that, without KLT truncation, the number of PCE coefficients would have been ∼ 1.3 • 10 6 , thus making the model training unfeasible.

Mean and standard deviation of the transfer functions are computed using [START_REF] Eldred | Recent advances in non-intrusive polynomial chaos and stochastic collocation methods for uncertainty analysis and design[END_REF] and [START_REF] Weng | Stochastic simulation and robust design optimization of integrated photonic filters[END_REF] and the absolute values are reported in Fig. 3(a) and (b) for the two ports using red solid and dashed lines, respectively. PDFs of the transmission at the two ports at the central wavelength, efficiently computed by sampling the model with the same samples used for the reference Monte Carlo analysis, are shown with red solid lines in Fig. 3(c) and(d). Results are in excellent agreement.

Lastly, we train the sparse PCE model described in Section II-C for the same set of 5 independent variables resulting from the truncation of the KLT. For the adaptive training, we set the expansion order to be in the range [START_REF] Bogaerts | Programmable photonic circuits[END_REF][START_REF] Wagner | Monte carlo simulation of electromagnetic scattering from two-dimensional random rough surfaces[END_REF], and we use only L = 100 samples, again drawn from a LHS design. The use of PCA with a truncation threshold of ε PCA = 10 -3 allows to vastly compress the 2000 wavelength samples (1000 for each of the two ports, bar and cross) to a mere ñ = 9 principal components. Individual sparse PCE models needs hence to be trained only for these principal components. It is worth noting that the UQLab toolbox only handles real data, therefore a separate model is trained for the real and imaginary part of the principal components, leading to a total of 18 models to be computed. The number of non-zero coefficients identified for these models by the adaptive algorithm ranges from 13 to 38. Results for the mean and standard deviation of the transmission are reported in Fig. 3 A similar analysis is performed for the seventh-order Mach-Zehnder filter, comprising 38 random variables. The probability density functions of the bar and cross transmissions at λ = 1521.8 nm are reported in Fig. 4(a) and (b), respectively. Blue bars refer to the results of the Monte Carlo analysis, red solid lines to the classical PCE computed after the application of the truncated KLT, and yellow dashed lines to the proposed sparse method in conjunction with PCA. By using the same truncation thresholds for the KLT and PCA as in the previous test case, the initial set of correlated parameters is reduced to d = 7 independent variables, whereas the output wavelength samples is compressed once again to only ñ = 9 principal components. An order of p = 6 is again used for the classical PCE, leading to K = 1716 coefficients. Also in this case, the training of the classical PCE model would have been unfeasible without the preliminary KLT truncation. Such a reduced model is instead trained with L = 3K = 5184 samples. For the adaptive sparse scheme, the same (maximum) order is used, and the models of the principal components is trained with only L = 200 samples. In both cases, training samples are generated from a LHS design. The number of non-zero coefficients identified ranges from 15 to 120. As can be seen in Fig. 4, similarly to the fifth-order filter, the transmission at the bar port is only minimally impacted by the variability, while the rejection at the cross port varies across a wide range. Here the most probable rejection value is -17 dB, compared to a value of -22 dB for the ideal case without parameter uncertainty.

To further demonstrate the accuracy of the proposed method, Fig. 5 (a) and (b) report the probability density functions of the 3-dB bandwidth for the fifth-and seventhorder filters, respectively. The bandwidth is derived as a postprocessing step starting from the full transfer functions described above, computed with a classical (solid red lines) and sparse PCE (yellow dashed lines) models after the application of the truncated KLT, and using PCA in the latter case. The bandwidth calculated from the Monte Carlo samples is used as reference (blue bars). Despite the fact that the analysis is conducted on a derived synthetic quantity (the bandwidth), the proposed method still retains its high accuracy, with an almost perfect agreement between the three techniques. It should also be noted that the two probability density functions are centered at different values because the designs of the two filters targets two different nominal bandwidths (39 GHz for the fifth-order filter and 29 GHz for the seventh-order one). Both density functions are symmetric, indicating an equal probability for the bandwidth to become either larger or narrower under the considered parameter uncertainty.

RMSE

Fifth-order filter (d = 28)

Seventh-order filter (d = 38) average maximum average maximum L = 25 2.0999×10 -2 3.0723×10 -2 2.6623×10 -2 4.7079×10 -2 L = 50 6.0095×10 -3 1.0539×10 -2 2.7191×10 -2 5.1511×10 -2 L = 100 3.7550×10 -3 5.8772×10 -3 9.1346×10 -3 1.5714×10 -2 L = 200 3.8724×10 -3 5.7735×10 -3 5.6308×10 -3 9.6840×10 -3 L = 400 2.7967×10 -3 4.1644×10 -3 3.7206×10 -3 5.8509×10 -3

D. Efficiency and Accuracy

Table I summarizes the efficiency of the discussed methods in terms of number of required samples and computational time for the analysis of the full transfer functions of the two test cases. All simulations were performed on a Lenovo ThinkPad X1 Yoga laptop with an Intel(R) Core(TM) i7-7500U processor, CPU running at 2.7 GHz, and 16 GB of RAM. It is important to mention that most of the computational time required by the PCE methods (about 98% for the classical method, and 65% for the proposed sparse method) is due to the simulation of the training samples.

As already noted, for both the reported examples, the application of the classical PCE would have been unfeasible without the preliminary KLT reduction, since it would have required a number of training samples in the order of several millions. It is also important to mention that this is not necessarily the case for the sparse method, thanks to its capability of identifying only the non-negligible terms. However, its computational efficiency also benefits from this dimensionality reduction, since the algorithms needs to search through a largely reduced sets of terms. Because of this, the performance of the method is expected to reduce if a large number of independent or weakly correlated parameters is included, or if the amount of variability is increased up to the point of requiring a model of prohibitive complexity. The advocated method achieves a further speed-up of about one order of magnitude w.r.t. the classical PCE method (8.9× and 18.8× for the first and the second test case, respectively), and an overall speed-up of 64.2× and 33.3× against Monte Carlo.

Table II investigates instead the impact of the number of training samples on the achieved accuracy. The table reports the average and maximum, over frequency and port, of the root-mean-square error (RMSE) computed w.r.t. the 10000 MC samples. It is observed that the RMSE steadily decreases by increasing the training set size L. In particular, the maximum RMSE drops below 10 -2 for L ≥ 100 in the first test case (fifth-order filter), and for L ≥ 200 in the second test case (seventh-order filter). The improvement becomes marginal beyond those values. A generally good practice is to select a number of training samples proportional to the number of uncertain parameters, e.g., L > 3 • d or L > 10 • d, where d is the original number of uncertain parameters, and d is the number of reduced variables after KLT.

IV. CONCLUSIONS

In this paper, we have demonstrated an efficient and accurate method to analyze the performance variability induced by fabrication uncertainty in photonic circuits with a large number of Gaussian correlated design parameters and multiple outputs of interest. The proposed scheme combines a sparse PCE approach with KLT to both remove correlation in the parameter set and reduce its dimensionality, alleviating the computational cost of the PCE model. Similarly, the use of PCA compression allows reducing the number of output quantities to be modeled, thereby further improving efficiency.

As an example of the method potentialities, we have shown the stochastic analysis of a cascaded Mach-Zehnder filter with up to 38 design parameters with a complex correlation scenario. The filter transfer function at different wavelengths represented the multiple output of interest for the analysis. Results were in excellent agreement with the reference Monte Carlo simulation, while the computational time was reduced by more than 30 times.

With the widespread of photonic devices controlled by large sets of parameters and with multiple objective requirements, this method represents a promising and valuable approach for the stochastic analysis, yield estimation, and optimization of their performance.

Fig. 1 .

 1 Fig. 1. Schematic of the fifth-order unbalanced Mach-Zehnder filter.

Fig. 2 .

 2 Fig. 2. Transfer function of the fifth-order (bold balck lines) and seventhorder (bold blue lines) Mach-Zehnder filter at the bar and the cross ports for the nominal design

Fig. 3 .

 3 Fig. 3. Stochastic analysis of the transfer function for the fifth-order Mach-Zehnder filter (28 random variables). Magnitude of the transfer function at the bar (a) and cross (c) port is considered. Solid gray lines represents a subset of the Monte Carlo samples. Other solid lines show the mean computed with Monte Carlo (blue), classical PCE (red), and sparse PCE with PCA (yellow). Dashed lines are the standard deviations computed with the same three methods. The probability density functions computed at λ = 1521.8 nm, dashed black line in (a) and (c), are reported as well for the bar port (b) and the cross port (d). The reference histogram (blue bars) is obtained from the Monte Carlo samples while solid red and dashed yellow lines report the results of classical PCE and sparse PCE with PCA, respectively.

Fig. 4 .

 4 Fig. 4. Stochastic analysis of the seventh-order Mach-Zehnder filter (38 random variables). Probability density functions of the transmission at λ = 1521.8 nm for the bar port (a) and cross port (b) are computed with Monte Carlo (blue bars), classical PCE (red solid lines) and the proposed sparse PCE with PCA (yellow dashed lines).

  (a) and (b) for the two ports using yellow solid and dashed lines, respectively, while probability density functions at the central wavelength are shown in Fig. 3(c) and (d) with yellow dashed lines. Results are once again in excellent agreement with the Monte Carlo analysis. However, compared to the classical PCE, a much smaller number of training samples is needed.

Fig. 5 .

 5 Fig. 5. Probability density function of the 3-dB bandwidth of the cross port at λ = 1521.8 nm. (a) Fifth-order filter and (b) seventh-order filter. The reference histogram is computed with Monte Carlo. Results for classical and sparse PCE are reported with solid red and dashed yellow lines, respectively.

TABLE I COMPARISON

 I OF THE DIFFERENT METHODS USED FOR THE ANALYSIS OF THE MACH-ZEHNDER FILTERS WITH 28 AND 38 VARIABLES.

		Number of samples Computational time
		d = 28	d = 38 d = 28	d = 38
	Monte Carlo	10000	10000	847 s	899 s
	Classical PCE	4 mln	21 mln n/a	n/a
	Classical PCE + KLT	1386	5148	118 s	507 s
	Sparse PCE + KLT + PCA 100	200	13.2 s	27 s
		TABLE II
	ACCURACY OF THE PROPOSED METHOD AS A FUNCTION OF THE NUMBER
		OF TRAINING SAMPLES.