Michel Labour

Tagging the interactive quality of Learning Objects

Keywords: multimedia, interactivity, CALL, metadata, tagging, document

published or not. The documents may come

Tagging the interactive quality of Learning Objects in CALL: the role of metadata in assessing multimedia for educational purposes

Introduction

One of the stumbling blocks of the New Technologies of Information and Communication is their apparent over-complexity that dissuades learners and teachers from having a real say in what has been, and is being, developed for them. This arises from the premise that a Tutoring System should take into account different ways of learning and helping learners to choose the best way to attain their objectives. This premise rests on the assumption that users understand the system enough to make appropriate choices. However, if there is no appropriate conceptual and practical access to Learning Objects (LOs 1 ) it is nigh impossible for users to utilise and to offer constructive feedback of their needs and wants in optimising the Tutoring System.

As the number of (learning) objects continues to grow exponentially and especially as our needs for learning expand equally dramatically, the lack of information or metadata about objects has produced a critical and fundamental constraint on our ability to discover, manage and use (learning) objects (Learning Technology Standards Committee. Learning Object Metadata. Working Draft Document 3, 1999-11-27).

In this context, Littlejohn (1998, page 205) speaks of the necessity of closely examining pedagogic materials so that learners and teachers can see "inside" them to take more control of their design and use. Indeed, descriptions of digital LOs in language learning that do exist, invariably, tend to focus on the more "passive" elements of the product (how many hours of video, linguistic level, technical shortcomings, etc.). Little is said, for example, of the degree of interactivity that the product offers, and even less about the different types of interactivity available or how these features are made known to potential users.

Based on the work of an interdisciplinary team at the University of Valenciennes (France) we will attempt to clarify this issue by proposing a system of classification with five different types of interactivity in computer-supported Tutoring Systems. Without claiming that our classification is exhaustive, our aim is to propose a framework of analysis to facilitate dialogue between users and creators in Computer Assisted Language Learning (CALL). In so doing, this study falls into the ambit of work in progress of organisations such as the Instructional Management System project of EduCause and the Learning Technology Standards Committee (see References below).

In this paper, our classification limits itself to a descriptive analysis of different types of interactivity, rather than an examination of the impact of technology on learners. The main objective of our typology is to facilitate the identification of various interactive elements by using metadata in a way that learners can make more informed decisions about the learning content and the way to learn it.

To do this, we will propose five broad categories of interactivity. We will start by briefly looking at the issue of tagging, followed by an interpersonal theory of learning which underlies our system of classification. After this, a second category of classification is presented involving ways that learner output is evaluated as evidence of cognitive interactivity. Then follows a three-part analysis using a scenistic approach to classify interactive features provided by the scenario. To illustrate the practical use of this classification three examples of multimedia products will be given. Finally, some implications for designing and developing CALL are proposed.

Tagging

The dilemma that tagging 2 LOs present is twofold. First, where does one draw the line regarding the level of detail that needs to be indicated regarding LOs i.e. a narrow, detailed or general, open-ended description? If the functional level of semantic content tagging is not appropriate, the degree of reusability of LOs could be severely limited. Second, should the vocabulary of these tags be fixed or essentially open-ended? On the one hand, a fixed vocabulary increases its interoperability. But on the other hand, too dense and rigid a vocabulary may limit its accessibility in terms of the different needs of authors, designers and users.

In attempting to answer these questions, we will look at tagging for "educational3 " use, notably that which concerns cd not limit itself to the flow of data between the user and LOs. It also includes the flow of data between the user and people-resources via the Tutoring System. In so doing, we hope to demonstrate the necessity of having a fixed vocabulary in that it offers a common reference point to facilitate dialogue between users and experts.

An interpersonal theory of learning

Based on his experimental work, Vygotsky (1934Vygotsky ( /1962, pages 101-104,116) , pages 101-104,116) showed that learner interaction with active external pedagogic agents, like teachers, can directly shape and transform the mental development of learners. For Vygotsky (1934/1962, page 8) effective learning necessarily involves mediators like teachers and various socially constructed (conceptual and mechanical) tools. In this sense the role of teachers and mentors is to provide a guiding hand at certain "critical" moments in the learning process (cf. the "zone of proximal development4 "). One of the key elements in this process is the way language, as a pivotal conceptual tool, is learnt by the individual through interpersonal experiences (Vygotsky 1934(Vygotsky /1962, page 137), page 137). This exchange implies some form of interactivity in such a way that learners need to adapt to the teaching/learning resources, just as these resources (teacher, a pedagogic tool) need to take into account the individual's learning process, notably the person's "zone of proximal development". In this way both learners and learning sources (real or virtual) act on each other. While the classical vygotskyean approach tends to emphasises more "vertical collaboration" between mentor and learner, a more recent approach includes "horizontal collaboration" in language learning via peer learning/teaching (Legutke and Thomas 1991, pages 276-279) without negating the place of the designated mentor. This theory of learning, notably the need for learner accompaniment, implies different types of interchange and collaboration to facilitate the learning process. These various interchanges are not mutually exclusive and can be complementary facets of the pedagogic environment.

Type 0: No help is given to the user (presumably because the target user is not expected to "need" it) Type 1: Built-in help -The designer anticipates the user's questions and integrates the solutions within the product itself based on interactions between the user and the machine. This help mode represents contact with a tutor in absentia. It can take the form of preestablished glossaries and a dictionary. In this way, the learner is limited to the product designer's choice of questions and solutions. Normally, it is up to users to find out how and when to best use these "help" features. In her book on CALL, Demaizière (1986, pages 503-507) provides practical guidelines in how to optimise this function.

Type 2: Peer contact -This can facilitate learning by (horizontal) collaborative and co-operative activities in which the externalising of ideas and feelings helps learners to discover alternative ways of satisfying their learning needs. For instance learning a language via email partnerships where learners with different native languages: exchange information, suggest corrections to each other's target language and ask for clarification about various elements of the target culture (Labour 1999). One advantage of this system is the possibility of selfmanagement, especially among more autonomous learners.

Type 3: Reactive contact -It is up to the learner to take the initiative to contact the tutor (learner-led vertical collaboration). Access to the tutor can be on-line, by email, in a chat room, or off-line. The learner must not only be aware of having a question to ask, but must also want to make the effort to find a solution with the help of someone with greater expertise. However, from a language learning point of view, a certain level of tolerance of ambiguity and incompleteness of knowledge is sometimes needed (Naiman et al 1978). This implies that learners need to know when a question is conducive to achieving their objectives and when it is not. Not to do this, may lead learners either not to ask any questions or to ask so many questions that they become side tracked from their initial learning objectives.

Type 4: Proactive contact -This type of interchange already exists in the traditional audio language laboratory where the teacher can take the initiative to intervene at the appropriate moment while the student is at work (tutor-led vertical collaboration). To do this, however, the tutor must first observe ("listen in" to) the learner at work and then step in at the right moment to encourage or to suggest corrections. While the approach may reassure some learners that they have not been left to their own devices, others may see this proactive intervention as tutor intrusion.

Having completed this first category of classification (see figure 3 below) it would now be useful to look at how one can assess the role of CALL from a more "cognitive" perspective.

Feedback on learner output

One of the aims of evaluating learner output is to assess the effectiveness of certain types of probable cognitive interactivity. In this sense, cognitive interactivity involves mental activity to do with data acquisition, storage, retrieval and use in a given context with which it interacts, as well as (reflexive) interactions with itself. In terms of learning, this implies that, on the one hand, learners can be encouraged to monitor the efficiency of their mental processes in achieving set goals (result-oriented metacognition). On the other hand, learners can be helped to be aware of their cognitive "tendencies" (e.g. use of memory, problem solving tactics) of processing new and old data (process-oriented metacognition). In light of this, a classification is proposed in evaluating output as a probable sign of significant cognitive interactivity.

Type 0: Read/listen -The learner receives the data but there is no associated activity to ensure that data have been successfully processed.

Type 1: Binary feedback -Each item of learner output (e.g. filling in blanks) is evaluated by a right or wrong response. No other comment or help is given to the learner.

Type 2: Template -This involves a more global process of comparing the learner's answer with that of the tutor's. For example learners can record their voices and compare it with a model (template). This can be done by voice recognition functions or ossilograms (to compare the learner's voice with that of the model) for oral expression, or by making available model answers to written work. Different variations of the template can be proposed in terms of the learner's needs.

Type 3: Record keeping -A global score is given for work done. This allows learners to have an idea of their rate of progress. This function can also include record keeping of work done, and that which needs to be done, i.e. an (automatic or semi-automatic) electronic logbook.

Type 4: Error correction -This involves giving learners a personalised and relatively detailed response to their errors. Each error has a different type of response to make the learner aware of why the mistake was probably made. This type of feedback tends to encourage process-oriented metareflection (see Demaizière 1986, pages 478-503).

Figure 1 below presents interactivity linked to learner accompaniment and feedback on output.

.

Role of the scenario

In an extensive study into designing effective multimedia texts, Leleu-Merviel (1996) proposes a "scenistic" analysis concerning the interactive quality of the scenario in determining what information is given to the user and the way this is done. In this context, a "scenario" is defined as a deliberately thought-out text that has: a story line, technical details, the author's intention, a target audience, and a message. From this point of view, the concept of "interactivity" implies a scenario which anticipates the various responses of the user that are oriented in a particular direction by a "scenational schema", or a series of pre-planned events/states linked together by different paths for the user to take. In this way with each interactive session the user can follow a different path. Thus when evaluating the "interactive" quality of CALL it is important to closely examine the choices the Tutoring System gives to the learner. In this context, Laubin, Durand and Leleu-Merviel (1999) propose a classification based along three axes:

• the interactive nature of the scenario, • the generative nature of the data, • the evolvibility of information units.

Before examining these categories more closely, it is important to underline that a product is never situated purely on one level. It may have different predominant features, but depending on the needs of the learning context other features may exist that can be classified into different categories. This, as we shall see below, is the case of three language learning cd-roms chosen for evaluation. We have seen how the concepts of interpersonal and significant cognitive interactivity can be operationalised, the challenge, however, is not simply in making these resources available in a haphazard way, but as part of a coherent whole, that is to say when a learning resource is most likely to be useful to the learner. It may be that different activities are of high technical excellence but lack the necessary overall interconnections to facilitate appropriate learning. Knowing that a learning tool has different pedagogic features, like voice recognition or a wave form oscillogram does not tell us how they are linked to other learning activities. Not to look at the overall coherence of the system, may lead to cognitive overload or to other forms of dysfunctions such as learners loosing sight of their priorities. One way that these learning resources can be co-ordinated to facilitate the learning process is by a carefully designed interactive scenario. In the discussion that follows, the role of the scenario will be examined in some detail.

6. The interactive nature of the scenario Type 0: Linear non-interruptible / audio-vision (e.g. cinema) -The system works independently of the reactions of the user (sit back and enjoy). While there is no direct interactivity with the user, there may be concomitant cognitive, emotional or social interaction.

Type 1: Linear interruptible / reading (e.g. video recorder, book) -Users can stop the text, skip sections and carry on as they please, or sometimes go back to reread the text. Essentially, however, the text remains linear and users have limited control in the way the content is presented.

Type 2: Menu-based navigation (e.g. arborescent database) -The navigation is done by a series of limited choices offered by a predetermined menu. The arborescent structure of the navigation offers little possibility for undoing the choice (unidirectional navigation), except if a previous choice has been "copied and pasted" ("cloned") in a successive phase.

Type 3: Multiple-criteria exploration (e.g. train time table) -Multi-criteria requests can be made, as at each step the choice is multidimensional. Users enter different criteria and then select what corresponds to their needs. In reality, it can be difficult to establish if a product is mainly type 2 or 3 as it may take too long (40-60 hours) to look at everything in the product. So one has to assess a sample of features.

Type 4: Free exploration (e.g. networks like the World Wide Web) -The exploration is less predictable. The user is not limited to predetermined paths. From a given entry point, networks of interconnected links allow the learner to create individualised paths that may not have been foreseen by the creator of the system. One can go back and forth on different paths that include segments of other greater routes. One can therefore jump in and out of different categories (help page, index page, news pages, etc.) unlike a strict arborescent route (type 2).

Type 5: Virtual visit (e.g. virtual/3-dimensional reality) -Users are fully immersed in a near continuous interactive universe, and choices are made by trial and error as users work their way through the system. The whole system depends on the continued action of users who act as if they are part of the same space as that of the electronic universe. In this type of 3-D space there is direct manipulation of objects without having to learn external metaphors. The quality is judged not by the richness of the "scenario" (different messages that can be sent and received) but by the quality of the "events" (data flow).

The generative nature of the data

The preparation of a "document" (i.e. an organised structure of content) involves the judicious choice of data and organising them in a given structure. It is the type of data and the way they have been organised that determines the different "scenational" modes.

Type 0: Fixed data -Choices and solutions are rigidly predetermined. There is no possibility of modifying data other than the way they were initially presented. Every time the application is used the user is exposed to the same content in exactly the same way.

Type 1: Basic automatic start-up generativity -There is an automatic, machine-generated selection of parameters at the start of the work session, independent of the action of the user. This feature adds a factor of initial unpredictability to the learning task. There is no way for the user to know on what basis the selection of the new data is done as the controlling variables are ensconced in the programming codes.

Type 2: Basic user-driven start up generativity -Users decide at the start of the work session which features and parameters they wish to activate (e.g. level of difficulty). One has access to some of the variables used, this makes evaluation of the learners' performance clearer than Type 1.

Type 3: Basic automatic dynamic generativity during use -The system automatically modifies the variables at the start and during the work session without the user being able to directly control the process. This can add an element of unpredictability (and realism) to the learning activities. That said, it also makes evaluation of progress more complex as it becomes much more difficult to quantify the interaction of the user with the tool.

Type 4: Basic user-driven dynamic generativity during use -The user can modify the parameters at the start and during the work session. The modification depends directly on the user's choice. The interaction is easier to measure than that of type 3.

The evolvibility of information units

It is often the case that data do not have meaning in themselves, but depend on their semantic links with other data. It can sometimes be necessary to unclip and reclip (regroup) different data together to achieve a given objective.

Type 0: Totally integrated fusion -Information units are completely integrated/indivisible both to its context and to the fact that they are one whole, with no sub-component parts, for example a traditional-format photograph.

Type 1: Indivisible & context-bound -Information units can be organised in various ways as long as they remains in their initial-state and context. One can reassemble the information units (e.g. displacing the different parts of a face in a specially formatted photograph) but one cannot break them down any further to create new (smaller, bigger, different) information units.

Type 2: Divisible & context-bound -Information units can be broken down into sub-component parts and reorganised in a different way but only within the initial context. For example, the different elements of a web page can be selected and copied onto another web page.

Type 3: Indivisible & transportable -Information units cannot be broken down any further but can be exported to other contexts that are different to the initial one, for example from a web page to a word processor page.

Type 4: Divisible & transportable -Information units can be broken down into their component parts and transported from their original context to a different context.

The role of the scenario in encouraging interaction between the user and the learning content can be represented on a graph in the following way (figure 2).

Some examples

In order to illustrate the use of this type of classification in CALL it would be useful to look at three language learning multimedia products currently used in learning English as a Foreign Language. The first cd-rom, The Adventures of KC in Londonc revolves around the theft of the Crown Jewels from the Tower of London. Learning activities are based on a variety of grammar, vocabulary, comprehension, and listening skills games, with an Avatar, called KC, who tries to get the jewels back from the Alphabet Gang. The second cd-rom, English Plus (version 2.0) contains 9 sections: reading, writing, language, vocabulary, web links, exploration of the cd-rom, and a virtual reality "adventure" game. It is this last feature that we will be examining. Finally we will look at, Top Score (International version) prepares students for the American Test Of English for International Communication (TOEIC) and one of its activities includes a 3-D Maze (vocabulary acquisition activity) which will be tagged according to our typology in Figure 3.

This type of description can be used to compare certain features of related multimedia products to help users decide which to use, or at least to be forewarned of their limits. By the same token, the description can serve as a common reference point when consulting users about this type of tutoring. Armed with these descriptors, the researcher can also set about evaluating the impact and the efficiency of a given Tutoring System. This could be done by conduc ting experiments in varying the different types of "interactivity". For example, one could compare learning outcomes and learner satisfaction via a Tutoring System composed of Scenation Type 5 5 + Generative Data Type 2 6 + Information Units Type 3 7 + Feedback on Output Type 3 8 + Learner Accompaniment Type 1 9 to those using a similar configuration but where one of the variables is changed, for example Learner Accompaniment Type 3 10 (rather than Type 1). These type of studies could help clarify key issues about what constitutes appropriate pedagogic accompaniment in terms of pedagogic content and presentation in CALL.

Conclusion and implications

At the risk of stating the obvious, it must be made clear that interactivity and user-control in CALL are not ends in themselves. There are times when limited interactivity and control for a given task within a given period of time may be more useful to the learner than a system that is too open ended to satisfy specific learning objectives. This begs the question, however, of knowing when to offer learners different types of interactivity. Hence the need for more research in this domain.

It is in this context that our classification attempts to provide common reference points, namely though agreed metadata, between system-designers and users concerning different types of interactivity and control that can be made available according to the demands of the learning objectives, pedagogic content, learner preferences and tutoring context. This interchange of perspective can foreshadow future innovations.

5 virtual visit 6 initial-user driven data 7 indivisible and context-bound information units 8 record keeping feedback 9 built-in help 10 reactive tutor-learner accompaniment With this in mind Tutoring Systems need to encourage learners to engage in motivating, exploratory activities that allow users to express their personality and in so doing develop more creative approaches to learning and using, in this case, English as a Foreign Language. This involves understanding the workings of the user's imagination and learning process via a scenational schema that creates semantic links between what the learner knows, wants to know, and what is not known. Ways of doing this include helping learners not only to acquire (declarative) knowledge but also how it can be used, managed and updated through a variety of tasks that encourage questioning, challenging, explaining, constructing, and tutoring.

From an educational perspective, CALL needs to more clearly define the limits of working by oneself and of interacting with others (peers, tutors) to accomplish a given task. The issue of the learning activity and context also includes practical factors like the resources available to learn, given that one of the characteristics of adult (language) learners is the need to learn "quickly and efficiently" (Labour 1998, page 155). This has direct implications in the way the Tutoring System helps learners choose appropriate pedagogic content and assesses their progress (e.g. Evaluation of Output Type 3 or Type 4?).

One of the challenges in designing a Tutoring System remains the judicious selection of appropriate elements in constructing the underlying model of how users interact with the system given that learners cannot be offered all they desire, if nothing else, for logistic reasons. The complexity of this task cannot be the exclusive province of any one academic discipline but that of interdisciplinary teams of researchers and designers. This is often easier said than done in that it demands a redefinition of the role of specialists from various related fields of research. For when it comes to modelling the user in a Tutoring System one faces the challenge of, for example, weighing up which of the different learners' needs and representations are to be included, or not, from the model. This raises complex issues, such as: How does one go about identifying the various needs and wants of learners in a virtual environment? Are certain needs and wants more of a hindrance than a help in achieving certain set objectives (see Labour 1998)? If so, how can one help the learner develop compensatory strategies? Is it to be done by the Tutoring System (Generative Data Type 3) or by learners themselves (Generative Data Type 4) through, for example, self-evaluating questionnaires?

Similarly, the concept of learner needs and goals are themselves problematic. Even when one has identified them, how can they be efficiently updated, given that learners' needs and objectives change with time. Does the system have an ongoing assessment of learners needs and goals linked to appropriate learning activities and tutoring styles? This implies looking more closely at factors that enhance pedagogic coherence of the teaching/learning experience, as well as at the density and timeliness of learning and tutoring activities.

In short, current research in the domain needs to broaden its scope of enquiry in examining how users can more easily identify, and consequently use, the interactive elements of CALL so to better control what they learn and how their learning is done. This control should include letting learners determine the degree of unpredictability, realism and complexity of the (virtual) language learning context. In this way, CALL can offer greater possibilities for learning and using a language in a way that makes the guided learning context more effective and authentic for learners and tutors alike.

Figure 1 .

 1 Figure 1. Feedback on output and learner accompaniment

Figure 2 .Figure 3 .

 23 Figure 2. The degree of interactivity in a virtual environment's scenario

This is one of the items cited by the Learning Technology Standards Committee along with other tags such as: LifeCycle, MetaMetaData, Annotation, etc.

Vygotsky (1934Vygotsky ( /1962 pages 14,59,85,94,104) pages 14,59,85,94,104) emphasises that pedagogic interaction is only effective if the learner is cognitively ready to accommodate new data. In terms of cognitive development, the "zone of proximal development" represents an area of possible growth in a give n time span. When learners are in this zone, the tutor can, with careful accompaniment, help them overcome their difficulties. This accompaniment can include an approach that, for example, progressively up-scales (scaffolding) the learner's role within a task until the person overcomes the obstacles unassisted. Once learners accomplish the task without the guidance of someone with greater expertise, it can then be concluded that significant learning has taken place.


Antoine Bambade
email: antoine.bambade@inria.fr

Sarah El-Kazdadi

Adrien Taylor

Justin Carpentier

ProxQP: Yet another Quadratic Programming Solver for Robotics and beyond

Quadratic programming (QP) has become a core modelling component in the modern engineering toolkit. This is particularly true for simulation, planning and control in robotics. Yet, modern numerical solvers have not reached the level of efficiency and reliability required in practical applications where speed, robustness, and accuracy are all necessary. In this work, we introduce a few variations of the well-established augmented Lagrangian method, specifically for solving QPs, which include heuristics for improving practical numerical performances. Those variants are embedded within an open-source software which includes an efficient C++ implementation, a modular API, as well as best-performing heuristics for our test-bed. Relying on this framework, we present a benchmark studying the practical performances of modern optimization solvers for convex QPs on generic and complex problems of the literature as well as on common robotic scenarios. This benchmark notably highlights that this approach outperforms modern solvers in terms of efficiency, accuracy and robustness for small to medium-sized problems, while remaining competitive for higher dimensions.

I. INTRODUCTION

Over the past decades, optimization has become an essential component of robotics and a key enabler to simplify and systematize the programming of complex robot movements. Nowadays, many robotic problems -ranging from simulation, control, planning and estimation -are framed as optimization problems. An important class of such optimization problems is that of convex quadratic problems, which allows dealing with, among others, friction-less unilateral contact modelling [START_REF] Redon | Gauss least constraints principle and rigid body simulations[END_REF], constrained forward dynamics [START_REF] Carpentier | Proximal and sparse resolution of constrained dynamic equations[END_REF], inverse kinematics and dynamics for task control [START_REF] Escande | Hierarchical quadratic programming: Fast online humanoid-robot motion generation[END_REF][START_REF] Herzog | Momentum control with hierarchical inverse dynamics on a torquecontrolled humanoid[END_REF][START_REF] Kuindersma | Optimization-based locomotion planning, estimation, and control design for the atlas humanoid robot[END_REF], and legged locomotion [START_REF] Wieber | Trajectory free linear model predictive control for stable walking in the presence of strong perturbations[END_REF], to name a few. QPs are also commonly used also as subroutines for solving more complex problems as for instance in the context of constrained optimal control problems [START_REF] Leineweber | An efficient multiple shooting based reduced SQP strategy for large-scale dynamic process optimization. part 1: theoretical aspects[END_REF][START_REF] Houska | Acado toolkit-an open-source framework for automatic control and dynamic optimization[END_REF][START_REF] Tassa | Controllimited differential dynamic programming[END_REF].

Formally, a QP corresponds to the minimization of a convex quadratic cost under some linear equality and inequality constraints. It is mathematically described as:

min x∈R d 1 2
x T Hx + g T x

s.t. Ax = b, Cx ≤ u, (QP)
where H ∈ R d×d is a symmetric positive semi-definite matrix (notation H ∈ S d + ), A ∈ R ne×d , C ∈ R ni×d , b ∈ R ne , and u ∈ R ni . d is the so-called problem dimension, while n e and n i are the numbers of equality and inequality constraints respectively. In many scenarios, QP instances have to be solved at very high-frequency (e.g., 1 kHz for inverse dynamics), under various levels of accuracy depending on the application, and potentially in relatively large dimension (e.g., model predictive control). Reliable, accurate, fast and numerically robust solvers are thus essential in modern robotic applications. Several methods and their associated solvers have been developed by the optimization community (e.g., qpOASES [START_REF] Ferreau | qpOASES: A parametric active-set algorithm for quadratic programming[END_REF], Quadprog [START_REF] Goldfarb | A numerically stable dual method for solving strictly convex quadratic programs[END_REF], OSQP [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF], Gurobi [START_REF]GUROBI optimizer reference manual[END_REF], Mosek [START_REF] Mosek | MOSEK optserver documentation[END_REF], BPMPD [START_REF] Mészáros | The bpmpd interior point solver for convex quadratic problems[END_REF], IPOPT [START_REF] Wächter | On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming[END_REF], OOQP [START_REF] Gertz | Object-oriented software for quadratic programming[END_REF], etc.). Interestingly, in the robotics community, Pandala et al. [START_REF] Pandala | qpSWIFT: A Real-Time Sparse Quadratic Program Solver for Robotic Applications[END_REF] have also proposed a dedicated solver based on off-the-shelf interior-point methods combined with adapted sparse routines, which has shown to operate efficiently for generating different gaits on quadrupedal robots.

To the best of our knowledge, there is currently no numerical QP solver meeting all the aforementioned requirements. In this work, we propose a new approach, based on the wellestablished augmented Lagrangian formalism [START_REF] Powell | A method for nonlinear constraints in minimization problems[END_REF][START_REF] Hestenes | Multiplier and gradient methods[END_REF][START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF], for solving generic QPs, and provide a corresponding efficient numerical C++ implementation. In this approach, we notably propose to combine the bounded constraint Lagrangian (BCL) globalization strategy with a (primal-dual) proximal method of multipliers applied to (QP). These contributions pave the way towards more advanced numerical methods for dealing with complex optimization problems in robotics, with the ambition of significantly reducing the computational burden, increase the numerical robustness of the solver while also lowering the need of manual tuning of the underlying hyperparameters. To validate our approach, we benchmark the new solver against against its state-of-the-art competitors on various problems ranging from randomly generated QPs to the hard Maros-Mészàros [START_REF] Maros | A repository of convex quadratic programming problems[END_REF] QPs of the optimization literature, which includes classical robotic problems.

The paper is organized as follows. In Section II, we review different commonly used methods for solving QPs. Section III defines the different notations, reviews the classic augmented Lagrangian techniques and introduces the proposed primaldual method for solving QPs. The core aspects of the method are detailed in Section IV. The complete algorithm is depicted in Section V. Finally, we benchmark the new solver against the state-of-the-art in Section VI.

II. RELATED WORK

Generic constrained convex QPs are commonly solved with iterative (or indirect) methods. These methods are traditionally divided into two main families: (i) active-set methods, and (ii) penalization methods. When the problem contains no inequality constraint, direct methods can also be used, as solving the corresponding QP can be done directly by solving the corresponding (linear) Karush-Kuhn-Tucker (KKT) system [START_REF] Nocedal | Numerical optimization[END_REF].

Iterative methods of both families typically aim at solving a cascade of simpler intermediary optimization problems whose solutions tend to those of the original problem.

1) Active-set methods: This family of methods, developed in the 50s, aims at determining the set of active constraints at an optimal point of (QP), see, e.g., [START_REF] Nocedal | Numerical optimization[END_REF]Section 16.5] for an introduction. Once this set of active constraints is determined, the constrained QP can be solved directly by considering a QP with equality constraints only, and whose solution matches that of the original QP. Popular active set-based convex QP solvers include the open-source qpOASES [START_REF] Ferreau | qpOASES: A parametric active-set algorithm for quadratic programming[END_REF], Quadprog [START_REF] Goldfarb | A numerically stable dual method for solving strictly convex quadratic programs[END_REF], and the QPA module in the open source software GALAHAD [START_REF] Gould | GALAHAD user documentation[END_REF].

On the negative side, active-set methods typically suffer from undesirable effects, such as "active-set cycling" [29, sections 13.5 and 16.5]. This kind of phenomenon might slow down the method when the QP at hand does not satisfy constraint qualification properties (such as the traditional "linear independence constraint qualification" (LICQ), see, e.g., [START_REF] Nocedal | Numerical optimization[END_REF]Definition 12.4]), that are sometimes hard to verify in practice and are often not met.

On the positive side, warm-start strategies can easily be incorporated within active-set methods. This is usually key for solving cascades of similar QPs, which is common in applications such as sequential quadratic programming (SQP) and model predictive control (MPC).

2) Penalization methods: This family of methods transform the original constrained problem into a sequence of problems with either no constraints or very simple ones (such as sign constraints on some variables). The new problems typically have objectives that consist in two terms: (i) the objective of the original problem, and (ii) a term for penalizing points not being feasible for the original problem. There are two very common types of penalization methods used in practice.

a) Interior-point methods: This family of penalization methods forces through a barrier function the sequence of intermediary optimization problems to have strictly feasible solutions with respect to the domain of the original problem; see, e.g., [START_REF] Nocedal | Numerical optimization[END_REF]Section 16.6] for an introduction. Primal-dual interior-point methods [START_REF] Mehrotra | On the implementation of a primal-dual interior point method[END_REF][START_REF] Wright | Primal-dual interior-point methods[END_REF] became popular in the 90s due to their good practical performances across a wide range of problems. Standard solvers using an interior-point method are commercial solvers Gurobi [START_REF]GUROBI optimizer reference manual[END_REF] and Mosek [START_REF] Mosek | MOSEK optserver documentation[END_REF], closedsource BPMPD [START_REF] Mészáros | The bpmpd interior point solver for convex quadratic problems[END_REF], open-source solver OOQP [START_REF] Gertz | Object-oriented software for quadratic programming[END_REF] and qpSWIFT [START_REF] Pandala | qpSWIFT: A Real-Time Sparse Quadratic Program Solver for Robotic Applications[END_REF].

Because of the homotopy-based structure of this family of methods, one of their main drawbacks is the difficulty to use them with warm-starting procedures, when solving a sequence of related QPs (e.g., within SQP or MPC settings).

b) Augmented Lagrangian-type methods: This family of methods is primarily based on the idea of Lagrangian relaxations with an additional quadratic penalization term (possibly piecewise) for encouraging feasibility of the iterate (see, e.g., [START_REF] Nocedal | Numerical optimization[END_REF]Section 17.3]). This kind of techniques emerged in the 70s through the works of Hestenes and Powell [START_REF] Powell | A method for nonlinear constraints in minimization problems[END_REF][START_REF] Hestenes | Multiplier and gradient methods[END_REF] and then later with those of Rockafellar [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF] which tightly emphasized their linked with the so-called "proximal-point method". Indeed, augmented Lagrangian-type methods naturally arise by applying proximal-point methods on either the dual or saddlepoint formulations to (QP), thereby offering the advantage of converging under relatively weak assumptions. Augmented Lagrangian-type methods also have the advantageous property of being able to exploit warm-start procedures. However, they might exhibit slow convergence behaviors in practice. Standard solvers based on augmented Lagrangians for solving QPs include the recent OSQP [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF] (via an alternating direction method of multipliers) and QPALM [START_REF] Hermans | QPALM: A Proximal Augmented Lagrangian Method for Nonconvex Quadratic Programs[END_REF].

III. PROXIMAL METHOD OF MULTIPLIERS FOR QUADRATIC PROGRAMMING

This section introduces the main contribution of the paper: a generic algorithm to solve QP problems of the form of (QP). We start by recalling the optimality conditions related to (QP). Then, we review the notions of augmented Lagrangian and proximal-point algorithms in the simpler context of equality constrained QPs. These concepts are then adapted to properly handle QP problems composed of both equality and inequality constraints. We conclude this section with a few practical considerations for rendering the proposed approach numerically more stable and practically more efficient.

A. Optimality conditions

The Lagrangian L associated to (QP) is defined by:

L(x, y, z) := 1 2 x T Hx + g T x + y T (Ax -b) + z T (Cx -u), (1)
with x ∈ R d , y ∈ R ne , z ∈ R ni + .
For linearly constrained convex optimization problems such as (QP), strong duality holds and the associated KKT conditions are necessary and sufficient for ensuring a primal-dual point (x, y, z) to be optimal (see, e.g., [1, Section 5.2.3] and [39, Section 2, page 5] for more details). For (QP), the KKT system is given by the set of equations:

       Hx + g + A T y + C T z = 0, Ax -b = 0, Cx ≤ u, z ⊙ [Cx -u] = 0, (KKT)
where ⊙ denotes the Hadamard product (i.e., for two vectors

u, v ∈ R d , u ⊙ v ∈ R d is the vector whose ith entry is u i v i ).
In practice, we look for a triplet (x, y, z) satisfying these optimality conditions (KKT) up to a certain level of predefined accuracy ϵ abs > 0 (dependent of the application), leading us to the following natural absolute "stopping criterion":

   ∥Hx + g + A T y + C T z∥ ∞ ≤ ϵ abs , ∥Ax -b∥ ∞ ≤ ϵ abs , ∥[Cx -u] + ∥ ∞ ≤ ϵ abs . (2)
The ℓ ∞ norm is preferred to the ℓ 2 norm as it is independent of the problem dimensions. It is also common to consider relative convergence criteria for early-stopping, as absolute targets might not be reached due to numerical issues [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF][START_REF] Mosek | MOSEK optserver documentation[END_REF][START_REF]GUROBI optimizer reference manual[END_REF].

B. Equality-constrained quadratic programs

In this section, we provide a high-level overview on the proximal method of multipliers (PMM), which is at the heart of our approach. PMM is closely related to the probably even more famous augmented Lagrangian method (ALM), which we review first in the context of equality-constrained QPs (i.e., when n i = 0).

a) Augmented Lagrangian: The ALM [36, Section 4] relies on augmenting the standard Lagrangian L, defined in (1), with a squared ℓ 2 penalization of the linear constraints:

L A (x, y; µ e ) := L(x, y, 0) + 1 2µe ∥Ax -b∥ 2 2 , (3)
where µ e > 0 is a positive penalty parameter. The ALM alternates between the minimization of L A with respect to the primal variables x and a simple update rule for the dual variables y:

x k+1 = arg min x L A (x, y k ; µ e ),
y k+1 = y k + 1 µe (Ax k+1 -b). (4)
Using the expression of L A , one can obtain a more explicit expression for the iterations in the case of equalityconstrained QPs:

x k+1 = -H + 1 µe A T A -1 [g + A T (y k -µ e b)], y k+1 = y k + 1 µe (Ax k+1 -b). (5)
At this stage, there are two conflicting goals to balance in this iterative process. First, the smaller the value of µ e the faster the convergence: indeed, 1 µe can be seen as a step-size of a proximal point method applied on the dual of the QP, as (4) can also be equivalently written as:

y k+1 = arg max y min x L(x, y, 0) -µe 2 ∥y -y k ∥ 2 2 , (6)
(see [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF]Section 4] for details), and the x k+1 update might be chosen as

x k+1 ∈ arg min x L(x, y k+1 , 0) -µe 2 ∥y k+1 -y k ∥ 2 2 . (7)
Second, as µ e gets smaller, the conditioning of (H + 1 µe A T A) gets worse, thereby limiting the numerical applicability of the approach, particularly when the condition number of A, thus the one of A T A (as cond(A T A) = cond(A) 2 ), is already potentially large [START_REF] Nocedal | Numerical optimization[END_REF]Section 17.1]. Such cases are typical in robotics, for instance in the context of kinematic singularity or redundant constraints, just to name a few. b) Proximal method of multipliers: PMM is an alternative to ALM with an additional proximal term on the primal variables [36, Equation 1.9], following the scheme:

(x k+1 , y k+1 ) = arg minmax x,y L(x, y, 0) -µe 2 ∥y -y k ∥ 2 2 + ρ 2 ∥x -x k ∥ 2 2 , (8)
where ρ > 0 is an additional penalty parameter ( 1 ρ corresponds to a step-size for the primal proximal term). Explicit maximization with respect to the dual variable allows recovering a method similar in spirit with the ALM formulation (4):

x k+1 = arg min x Φ k ρ (x), y k+1 = y k + 1 µe (Ax k+1 -b), (9)
where

Φ k ρ (x) := L A (x, y k ; µ e ) + ρ 2 ∥x -x k ∥ 2
2 is often referred to as the proximal augmented Lagrangian (PAL) [START_REF] Hermans | QPALM: A Newton-type Proximal Augmented Lagrangian Method for Quadratic Programs[END_REF]. In practice, one can directly solve [START_REF] Escande | Hierarchical quadratic programming: Fast online humanoid-robot motion generation[END_REF] via its optimality conditions, encoded in the following linear system of equations:

H + ρI A T A -µ e I x k+1 y k+1 = ρx k -g b -µ e y k . (10)
It is worth mentioning that this linear system involves a matrix that is always nonsingular thanks to the two regularization terms ρ 2 ∥x -x k ∥ 2 2 and µe 2 ∥y -y k ∥ 2 2 . In other words, the problem ( 8) is always well-defined in the iterative process. As a comparison, one can notice that the linear system depicting optimality conditions for ALM might be singular, encoding the fact that an acceptable x k+1 might not be unique in [START_REF] Kazdadi | Equality Constrained Differential Dynamic Programming[END_REF]. As a consequence, intermediary computations involved in PMM are numerically more stable by construction. Yet, another advantageous feature of PMM is that it is guaranteed to converge to an optimal primal-dual pair (x k , y k ) → (x * , y * ) under relatively weak assumptions (existence of an optimal primal-dual pair with zero duality gap; see [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF]Theorem 7]). Under similar weak assumptions (existence of an optimal dual solution; see [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF]Theorem 4]), ALM is only guaranteed to converge on the dual variable y k → y * . It is also worth noticing that in the case of a QP involving only linear equality constraints, one iteration of either ALM or PMM can be cast as the solution of single linear system. We see in the next section that the situation changes markedly with inequalities.

C. Inequality constrained quadratic programs

In the presence of inequality constraints in (QP), a natural extension of (8) consists in iterating

(x k+1 , y k+1 , z k+1 ) = arg minmax x,y,z≥0 L(x, y, z) + ρ 2 ∥x -x k ∥ 2 2 -µe 2 ∥y -y k ∥ 2 2 -µi 2 ∥z -z k ∥ 2 2 . (11)
By explicit maximization in (y, z), one can also reach an equivalent formulation in terms of

   x k+1 = arg min x Φ k ρ (x), y k+1 = y k + 1 µe (Ax k+1 -b), z k+1 = [z k + 1 µi (Cx k+1 -u)] + (12)
where [.] + stands for the (componentwise) nonnegative part, and where the PAL Φ k ρ (x) is now defined using the augmented Lagrangian formulation for the problem involving inequalities [36, Eqs 1.4-1.5]:

L A (x, y,z; µ e , µ i ) := L(x, y, 0) + 1 2µe ∥Ax -b∥ 2 2 + 1 2µi [Cx -u + µ i z] + ∥ 2 2 -∥µ i z∥ 2 2 , (13)
Φ k ρ (x) := L A (x, y k , z k ; µ e , µ i ) + ρ 2 ∥x -x k ∥ 2 2 . ( 14
)
One can note that this PAL is a piecewise quadratic function.

When only equality constraints are involved in the QP, the PAL is a simple quadratic function, and the solution of the intermediary subproblems are obtained by solving a single linear system. The situation becomes a bit more subtle in the context of inequality constraints. For this reason, the practical algorithm that is proposed and investigated below is based on an approximate version of PMM (first proposed in [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF]):

(x k+1 , y k+1 , z k+1 ) ≈ ϵ k arg minmax x,y,z≥0 L(x, y, z) + ρ 2 ∥x -x k ∥ 2 2 -µe 2 ∥y -y k ∥ 2 2 -µi 2 ∥z -z k ∥ 2 2 ( 15
)
where ≈ ϵ k stands for requiring (x k+1 , y k+1 , z k+1 ) to be an ϵ k -approximate solution to the intermediate saddle-point subproblem. We control the accuracy level of this approximation via the following condition:

∥r k (x k+1 , y k+1 , z k+1 )∥ ∞ ≤ ϵ k , (16)
where r k is some nonlinear operator gathering the KKT conditions for the saddle-point subproblem at iteration k. In others words ∥r k (x, y, z)∥ ∞ = 0 iff (x, y, z) is the solution to [START_REF] Ph | A primal-dual augmented Lagrangian[END_REF]. For more details on r k , see Section IV-B.

For going further, we must address two important remaining aspects: (i) choosing appropriate rules for setting ϵ k (large values for ϵ k will obviously not lead to convergence of this numerical method) and µ k e and µ k i (appropriate tuning of those step-size rules critically impacts the practical performance of the method); (ii) computing suitable triplets (x k+1 , y k+1 , z k+1 ) satisfying ( 16) in an efficient way.

In the next section, we review the BCL strategy (originating from [START_REF] Conn | A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds[END_REF]) for dealing with (i). This strategy is key in our practical implementation of the QP solver. The problem (ii) of finding suitable approximations for the proximal subproblems is handled in Section IV.

D. BCL globalization strategy

In this section, we explain our strategy for fixing the hyper-parameters of the solver (tolerance on subproblems ϵ k , step-sizes µ e and µ i ). We rely on BCL (see [START_REF] Conn | A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds[END_REF] and [START_REF] Nocedal | Numerical optimization[END_REF]Algorithm 17.4]) which has been proved to perform well in advanced optimization packages such as LANCELOT [START_REF] Conn | Lancelot a Fortran Package for Large-Scale Nonlinear Optimization (Release A)[END_REF] and also in robotics for solving constrained optimal control problems [START_REF] Plancher | Constrained unscented dynamic programming[END_REF][START_REF] Kazdadi | Equality Constrained Differential Dynamic Programming[END_REF][START_REF] Jallet | Implicit differential dynamic programming[END_REF]. For solving QPs, we propose to rely on the combination of BCL with the proximal method of multipliers from Section III-C.

The main idea underlying BCL consists in updating the dual variables y k and z k obtained from [START_REF] Hermans | QPALM: A Newton-type Proximal Augmented Lagrangian Method for Quadratic Programs[END_REF] only when the corresponding primal infeasibility (denoted by p k+1 hereafter) is small enough. More precisely, we use a second sequence of tolerances denoted by ϵ k ext (which we also tune within the BCL strategy) and update the dual variables only when p k+1 ≤ ϵ k ext , where p k+1 denotes the primal infeasibility as follows:

p k+1 := max(∥Ax k+1 -b∥ ∞ , ∥[Cx k+1 -u] + ∥ ∞ ). (17)
It remains to explain how the BCL strategy chooses appropriate values for the hyper-parameters ϵ k , ϵ k ext , µ i and µ e . As for the update of the dual variables, it proceeds in two stages:

• If p k+1 < ϵ k ext :
the primal infeasibility is good enough, we thus keep the constraint penalization parameters as is.

• Otherwise: the primal infeasibility is too large, we thus increase quadratic penalizations terms on the constraints for the subsequent proximal subproblems [START_REF] Hermans | QPALM: A Newton-type Proximal Augmented Lagrangian Method for Quadratic Programs[END_REF]. Concerning the accuracy parameters ϵ k and ϵ k ext , the update rules are more technical and the motivation underlying those choices is to ensure global convergence: an geometric-decay type update when primal infeasibility is good enough, and see [START_REF] Conn | A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds[END_REF]Lemma 4.1] for when the infeasibility is too large. The detailed strategy is summarized in Algorithm 1.

Remark 1 (Global convergence). The BCL strategy [4, Algorithm 1 and Algorithm 2] was originally developed using an augmented Lagrangian, as opposed to a proximal augmented Lagrangian in our case. As a consequence, some of the original convergence guarantees from [START_REF] Conn | A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds[END_REF] do not hold anymore. For this reason, we introduced a safeguard parameter k max ∈ N enforcing the algorithm to ultimately always accept the candidates multipliers from the PMM updates [START_REF] Hermans | QPALM: A Newton-type Proximal Augmented Lagrangian Method for Quadratic Programs[END_REF]. This simple trick allows to inherit some nice properties from PMM that include convergence under mild assumptions [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF]Theorem 4].

In the next section, we review the last missing piece in our approach: a method for approximately solving the intermediary proximal subproblems [START_REF] Hermans | QPALM: A Newton-type Proximal Augmented Lagrangian Method for Quadratic Programs[END_REF].

IV. SOLVING THE PROXIMAL SUBPROBLEMS

For obtaining an approximate solution to [START_REF] Hermans | QPALM: A Newton-type Proximal Augmented Lagrangian Method for Quadratic Programs[END_REF], we introduce a primal-dual merit function, originating and extending Gill and Robinson's primal-dual augmented Lagrangian (PDAL) [START_REF] Ph | A primal-dual augmented Lagrangian[END_REF] to the case of inequality constraints. We show that associated semi-smooth Newton steps involve linear systems, whose structure is well-conditioned and similar to the one presented for the equality constrained case [START_REF] Gertz | Object-oriented software for quadratic programming[END_REF].

In the next subsections, we first review the classic primal PMM (IV-A), before providing details on the proposed primaldual approach (IV-B).

A. A primal semi-smooth approach

One possible instantiation of ( 15)-( 16) consists in maximizing the PAL analytically w.r.t. the dual variables (y, z) and minimizing it approximately w.r.t. the primal variable x:

   x k+1 ≈ ϵ k arg min x Φ k ρ (x), y k+1 = y k + 1 µe (Ax k+1 -b), z k+1 = [z k + 1 µi (Cx k+1 -u)] + , (18)
where ≈ ϵ k now stands for x k+1 to be an ϵ k -approximate solution in the following sense:

∥∇ x Φ k ρ (x k+1 )∥ ∞ ≤ ϵ k . ( 19
)
The PAL function Φ k ρ is semi-smooth [START_REF] Mifflin | Semismooth and semiconvex functions in constrained optimization[END_REF]. Hence its unique minimum can be found in finite time using a semi-smooth Newton method with exact line-search (see, e.g., convergence proof in [40, Theorem 3] and algorithm in [15, Section IV.C]). Practically speaking, the semi-smooth Newton method is initialized at x(0) := x k and generates a sequence x(1) , x(2) , . . . via the update rule:

x(l+1) = x(l) + α * dx, (20)
where the step-size α * is computed via an exact line-search:

α * := arg min α≥0 Φ k ρ (x (l) + αdx), (21)
(note that α -→ Φ k ρ (x (l) + αdx) is a continuous piecewise quadratic function with a finite number of breaking points), and where dx is found by solving a linear system of equations:

  H + ρI A T C T I k (x (l) ) A -µ e I 0 C I k (x (l) ) 0 -µ i I     dx dy dz   =   -∇ x Φ k ρ (x (l) ) 0 0   , ( 22
) where I k (x (l) ) := {i ∈ [1, n i ] | C i
x(l) -u i + z k i µ i ≥ 0} refers to the active set of the current subproblem, and C I k (x (l) ) corresponds to a reduced version of the matrix C containing the active rows indexed by I k (x (l) ).

This iterative process is repeated until reaching the accuracy requirement [START_REF] Houska | Acado toolkit-an open-source framework for automatic control and dynamic optimization[END_REF]; that is, as soon as ∥∇Φ k ρ (x (l) )∥ ∞ ≤ ϵ k it outputs x k+1 ← x(l) as an approximate solution.

One drawback of such strategy comes from the right hand side term ∇ x Φ k ρ (x). Indeed, similarly to (5), ∇ x Φ k ρ (x) gathers quadratics involving the constraint matrices A and C, which ill-conditioning impact the whole accuracy in linear system [START_REF] Kuindersma | Optimization-based locomotion planning, estimation, and control design for the atlas humanoid robot[END_REF]. It can lead to numerical saturation effects with illconditioned problems when high precision is required. For this reason, we have decided to investigate other merit functions in order to solve better conditioned linear systems. It leads to our primal-dual approach.

B. A primal-dual approach

In this section, we present a primal-dual approach to solve [START_REF] Hermans | QPALM: A Newton-type Proximal Augmented Lagrangian Method for Quadratic Programs[END_REF]. We introduce a primal-dual merit function, originating and extending Gill and Robinson PDAL to the case of inequality constraints [START_REF] Ph | A primal-dual augmented Lagrangian[END_REF]. The associated semi-smooth Newton steps involve better-conditioned linear systems with a structure similar to the one presented in [START_REF] Gertz | Object-oriented software for quadratic programming[END_REF]. a) Optimality conditions: A triplet (x, y, z) is a solution to [START_REF] Ph | A primal-dual augmented Lagrangian[END_REF] if and only if [START_REF] Marchi | On a primal-dual Newton proximal method for convex quadratic programs[END_REF]:

   Hx + g + ρ(x -x k ) + A T y + C T z = 0, µ e y -(Ax -b + µ e y k ) = 0, µ i z -[Cx -u + µ i z k ] + = 0. (23)
Hence we define r k (x, y, z) as:

r k (x, y, z) :=   Hx + g + ρ(x -x k ) + A T y + C T z µ e y -(Ax -b + µ e y k ) µ i z -[Cx -u + µ i z k ] +   . (24)
b) Generalized primal-dual augmented Lagrangian: Reinjecting appropriately second and third equations of ( 23) in the first one, one can notice, that ( 23) is equivalent to [START_REF] Mehrotra | On the implementation of a primal-dual interior point method[END_REF]:

           Hx + g + ρ(x -x k ) + A T ( 1 µe (Ax -b) + y k ))+ C T [ 1 µi (Cx -u) + z k ] + + (Ax -b -µ e (y -y k )) +([Cx -u + µ i z k ] + -µ i z) = 0, µ e y -(Ax -b + µ e y k ) = 0, µ i z -[Cx -u + µ i z k ] + = 0. (25
) Condition ( 25) correspond to KKT conditions for the following primal-dual merit function M k µ,ρ :

M k µ,ρ (x, y, z) := Φ k ρ (x) + 1 2µ e ∥Ax -b -µ e (y -y k )∥ 2 2 + 1 2µ i ∥[Cx -u + µ i z k ] + -µ i z∥ 2 2 .
(26) This merit function is strictly convex and continuously differentiable. Hence, as for Φ k ρ in (IV-A), its unique minimum can be found in finite time using a semi-smooth Newton method with exact line-search [START_REF] Sun | On piecewise quadratic Newton and trust region problems[END_REF]Theorem 3]. Thanks to the equivalence between ( 23) and ( 25), r k can consequently be used as a suitable stopping criterion for measuring solution required inexactness:

∥r k (x k+1 , y k+1 , z k+1 )∥ ∞ ≤ ϵ k . ( 27
)
Remark 2. Gill and Robinson introduced a generalized PDAL function [START_REF] Ph | A primal-dual augmented Lagrangian[END_REF] (x, y) -→ G k ρ,µ (x, y) for first tackling equality constrained problems. In the presence of inequalities, this PDAL function can be framed in its equality constrained form introducing a slack variable s satisfying the new equality constraint:

Cx -u -s = 0. ( 28
)
The minimization of G k ρ,µ w.r.t. variables x, y, z and s commutes. Considering the problem structure and following ideas from [START_REF] Marchi | On a primal-dual Newton proximal method for convex quadratic programs[END_REF], it can be shown that s and z can be directly deduced as functions of x or z. Consequently, the ordering of variables used to minimize G k ρ,µ defines different merit functions which may inherit from useful features of G ρ,µ in order to build an iterative procedure for solving [START_REF] Hermans | QPALM: A Newton-type Proximal Augmented Lagrangian Method for Quadratic Programs[END_REF]. For example, minimizing G k ρ,µ w.r.t. the slack variable s and re-injecting its optimal value s (depending of x and z) in the PDAL allows obtaining the merit function used in QPDO solver [START_REF] Marchi | On a primal-dual Newton proximal method for convex quadratic programs[END_REF]. Another possibility, consists in minimizing G k ρ,µ w.r.t. z, y and then to s and to re-inject all optimal values found (which are functions of x). In the latter case one retrieves then the PAL merit function introduced first by Rockafellar [START_REF] Rockafellar | Augmented Lagrangians and Applications of the Proximal Point Algorithm in Convex Programming[END_REF]. In the sequel, in order to get better conditioned linear systems, we decide to follow the same strategy but to relax optimal y and z values found, which leads to the new merit function M k ρ,µ . c) Primal-dual Newton semi-smooth steps: A semismooth Newton step applied to M k µ,ρ , and initiated at (x (0) , ŷ(0) , ẑ(0) ) = (x k , y k , z k ) involves finding for l ≥ 0 dw := (dx, dy, dz) such that:

∇ 2 M k µ,ρ (x (l) , ŷ(l) , ẑ(l) )dw + ∇M k µ,ρ (x (l) , ŷ(l) , ẑ(l) ) = 0, (29)
where ∇ 2 M k µ,ρ (x (l) , ŷ(l) , ẑ(l) ) stands for an element of its generalized Hessian [35, section 23]. It reads equivalently:

                       H + ρI + 1 µe A T A + 1 µi ( Ĉ(l) ) T Ĉ(l) -A T -( Ĉ(l) ) T -A µ e I 0 -Ĉ(l) 0 µ i I     dx dy dz   =     -(∇ x Φ ρ (x (l) ) + A T ( 1 µe (Ax (l) -d) + y k -ŷ(l) ) +C T ([ 1 µi (C x(l) -u) + z k ] + -ẑ(l) ) Ax (l) -d -µ e (ŷ (l) -y k ) [C x(l) -u + µ e z k ] + -µ i ẑ(l)     , (30)
where I k (x) is the active set of the current subproblem at x:

I k (x) = {i ∈ [1, n i ]|Cx -u + µ i z k ≥ 0}, (31)
and Ĉ(l) is a short-hand for denoting the generalized Jacobian of [Cx -u + µ i z k ] + at x (l) :

Ĉ(l) i := C i if (Cx (l) -u + µ i z k ) i ≥ 0, 0 otherwise, ( 32
)
where C i denotes the ith row of C. Remarking that for inactive constraints, i.e., i / ∈ I k (x (l) ), we have:

dz i = -(ẑ (l) ) i , (33)
the linear system (30) can hence be equivalently formulated as:

                         H + ρI A T C T I k (x (l) ) A -µ e I 0 C I k (x (l) ) 0 -µ i I I k (x (l) )     dx dy dz I k (x (l) )   = -    H x(l) + g + ρ(x (l) -x k ) + A T ŷ(l) + C T I k (x (l) ) ẑ(l) I k (x (l) ) Ax (l) -d -µ e (ŷ (l) -y k ) ([C x(l) -u + µ i z k ] + -µ i ẑ(l) ) I k (x (l) )    , dz I c k (x (l) )) = -ẑ (l) I c k (x (l) )) , (34)
where I c k (x) is defined as the set of inactive constraints at x. Eq. ( 34) has the same structure as [START_REF] Gertz | Object-oriented software for quadratic programming[END_REF] in the equality constrained case. The right hand side does not contain any term involving square matrices. The main differences come from the active set introduced by M k ρ,µ , which filters the only constraints that need to be taken into account in [START_REF] Redon | Gauss least constraints principle and rigid body simulations[END_REF]. d) Primal-dual line-search procedure: Once a semismooth Newton step (dx, dy, dz) has been obtained, the exact line-search procedure consists in finding the unique α * such that:

α * = arg min α≥0 M k ρ,µ (x (l) + αdx, ŷ(l) + αdy, ẑ(l) + αdz). (35
) Similarly to [START_REF] Jallet | Constrained differential dynamic programming: A primal-dual augmented lagrangian approach[END_REF] the function α → M k ρ,µ (x (l) + αdx, ŷ(l) + αdy, ẑ(l) + αdz) is a continuous piecewise quadratic function with a finite number of breaking points. Hence, one can compute α * exactly. Finally, the primal-dual semi-smooth Newton method initiated at (x (0) , ŷ(0) , ẑ(0) ) = (x k , y k , z k ) generates a sequence (x l , ŷl , ẑl ) via the update rule:

x(l+1) = x(l) + α * dx, ŷ(l+1) = ŷ(l) + α * dy, ẑ(l+1) = ẑ(l) + α * dz.

V. DETAILED APPROACH

The complete ProxQP procedure is provided in Algorithm 1. It contains a few additional practical enhancements, which are now detailed.

Preconditioning. ProxQP contains a preconditioning strategy, enhancing the overall numerical stability and convergence of the optimization process. The preconditioner used in our current implementation is often referred to as the Ruiz equilibration [START_REF] Ruiz | A scaling algorithm to equilibrate both rows and columns norms in matrices[END_REF]; see, e.g., [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF]Algorithm 2].

Initialization. Motivated by the fact that equality constraints are always active at an optimal solution, we use the following initialization for primal and dual variables:

x init y init = H + ρI A T A µ e I -1 -g b and z init = 0, (36)
which corresponds to the (primal-dual) solution to the corresponding QP where the inequality constraints were removed.

For the same reason, the default initial penalization is larger for equality constraints than it is for inequality constraints.

The choice µ e = µ i /100 appears to perform well.

Cold restarts. To solve hard QPs to high precision, we use an additional cold restart strategy. The idea consists in resetting the penalization parameters µ e and µ i as soon as both primal and dual feasibility appears to stall while the current values for µ e and µ i are smaller than a certain threshold.

VI. RESULTS

In this section, we detail the software implementation of Algorithm 1 and then benchmark the practical performance of this new solver by comparing it to existing state-of-the-art solvers on various problems ranging from randomly generated problems to the harder Maros-Mészàros QPs [START_REF] Maros | A repository of convex quadratic programming problems[END_REF], which includes a few classic robotic problems.

Software implementation. Our solver, referred to as ProxQP, is implemented in C++ and is extensively rooted on the generic-purpose Eigen library [14] for linear algebra. The current implementation is tailored for dense matrix operations, Algorithm 1: ProxQP Inputs:

• initial states: x 0 , y 0 , z 0 ,

• initial parameters: ϵ 0 ext , ϵ 0 , ϵ abs , ρ, µ e , µ i > 0 • hyper-parameters: µ f < 1, α bcl ∈ (0, 1), β bcl ∈ (0, 1), k max ∈ N, µ i,min , µ e,min > 0. Initialization:

• preconditioning (see Section V)

• optional initialization (see Section V) of x 0 ,y 0 , z 0 . while Stopping criterion (2) not satisfied do Compute (x, ŷ, ẑ) satisfying ( 16) (ϵ k -approximation to proximal subproblem ( 15)) using Section IV-B;

x k+1 = x if p k+1 < ϵ k ext OR k ≥ k max then ϵ k+1 = ϵ k µ i , ϵ k+1 ext = ϵ k ext µ β bcl i y k+1 = ŷ, z k+1 = ẑ else µ i ← -max(µ i,min , µ f µ i ) µ e ← -max(µ e,min , µ f µ e ) ϵ k+1 = ϵ 0 µ i , ϵ k+1 ext = ϵ 0 ext µ α bcl i y k+1 = y k , z k+1 = z k end k ← k + 1
Apply cold restart if conditions are met end Output: A (x k , y k , z k ) satisfying the ϵ abs -approximation criterion (2) for problem (QP). and leverages recent CPU architectures providing advanced vectorization mechanisms. Our new solver if freely available at https://github.com/Simple-Robotics/proxsuite and comes with an easy-to-use interface inspired from OSQP [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF].

Additionally, we have developed a dedicated LDLT Cholesky factorization to explicitly account for the specific features of the proposed approach. In particular, this Cholesky factorization implements advanced update routines to efficiently accounting for the change of active sets when solving [START_REF] Kuindersma | Optimization-based locomotion planning, estimation, and control design for the atlas humanoid robot[END_REF] and [START_REF] Kuindersma | Optimization-based locomotion planning, estimation, and control design for the atlas humanoid robot[END_REF], while lowering the overall memory footprint to maximize the performances. The new LDLT Cholesky is freely available within the same repository https://github.com/Simple-Robotics/proxsuite and will be hopefully integrated in Eigen.

Contrary to other solvers such as OSQP, qpSWIFT or qpOASES, our current implementations only relies on dense linear algebra routines (we plan to deeply exploit the sparsity of sparse problems in future works). Overall, as highlighted by the subsequent results, ProxQP performs better than modern sparse solvers over relatively sparse problems ranging from small to medium size (up to d = 1000) while remaining competitive for larger dimensions.

Benchmark scenarios. We benchmark our implementation on different types of QPs: equality and inequality constrained QPs, degenerate QPs with only inequality constraints, non-strictly convex QPs with only inequalities, as well as hard QPs from the Maros-Mészàros dataset [START_REF] Maros | A repository of convex quadratic programming problems[END_REF] with different levels of sparsity and dimensions. We also benchmark on typical QPs encountered in inverse kinematic and dynamic robotic problems. The benchmarks are available at https://github.com/Bambade/proxqp benchmark with an easyto-use interface inspired from the one proposed by OSQP [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF].

Our tests were carried in the following conditions:

• The level of accuracy required for termination is set arbitrarily to ϵ abs = 10 -9 (see criterion in [START_REF] Carpentier | The pinocchio c++ library: A fast and flexible implementation of rigid body dynamics algorithms and their analytical derivatives[END_REF]) so that to show across all these experiments which solvers manage to be the most accurate, the fastest, and the most robust within a large variety of QPs. • A time limit of 1000 seconds has also been set as in OSQP API [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF] for benchmarks (if a problem is not solved by a solver at such precision its "solving" time is set to this value in order to draw test plots). • For all the tests, we used the same set of parameters for Algorithm 1: µ e = 10 -3 , µ i = 0.1, µ f = 0.1, ρ = 10 -6 , α bcl = 0.1, β bcl = 0.9, µ i,min = 10 -8 , µ e,min = 10 -9 . Finally, ϵ 0 and ϵ 0 ext follow the initialization procedure proposed in [START_REF] Nocedal | Numerical optimization[END_REF]Algorithm 17.4], the safeguard k max = 10 6 (not used in practice).

Remark 3 (Accuracy choice). Choosing a low ϵ abs has a few non-negligible advantages when comparing solvers. Indeed, the fact a solver might not reach every desired level of accuracy (within the available finite precision limits) is typically due to either (i) an algorithm that has a (very) slow convergence, or (ii) somehow "inappropriate" underlying subroutines (including linear algebra ones) with a limited working precision/stability range, worsening the effect of finite precision. Hence, such low accuracy level reveals which solvers provide on the same time (i) an algorithm whose capabilities allow to reach high accuracy and (ii) a set of underlying numerical routines allowing to reach high precisions in reasonable times.

Benchmark setup. Benchmarks are performed with a standard laptop equipped with a relatively old CPU (Core i5 -5300U -5th Generation @ 2,3 GHz processor) and, for some specific benchmarks, also with a more recent CPU (Intel(R) Core(TM) i7-4790 CPU @ 3.60GHz) typically equipping clusters, illustrating potential gains of using a more recent architecture.

A. Random problems

Concerning random problems, three different levels of sparsity (0.15, 0.5, and 1) were used for generating the random matrices H, A and C in the spirit of the benchmark API from OSQP [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF]. The dimensions of those problems were picked from d = 10 to d = 1000. For each set of parameters (sparsity levels and dimensions), we generated 5 problem instances with different random seeds and averaged the running time of each algorithm on 10 consecutive runs (to limit the impact of loading costs). The results are provided in Figure 1 using bar plots (including the median, the minimal and maximal execution timings). We also report base statistics in Table I. In particular, Figure 1 shows that even in a sparse configuration (we recall that our solver uses a dense back-end), our solver is around 4 to 5 times faster than OSQP (the second best solver from our test-bed) for examples of dimensions d ≈ 50 (which is representative of typical robotic applications). When dimension grows to d ≈ 1000, the speed-up reaches almost an order of magnitude.

The failure rates observed in Table I for MOSEK and GUROBI solvers come from the fact they are not able to reach the desired precision ϵ abs = 10 -9 : solutions proposed are not accurate enough. Consequently, it also impacts their geometric means.

B. Degenerate pure inequality-constrained problems

Figure 2 provides the results of our numerical experiments when generating pure inequality-constrained QPs where the matrix H is positive definite but for which LICQ conditions are no longer satisfied (by duplicating the constraints), a common type of degeneracy. We observe from Figure 2 Solving times for sparse degenerate pure inequality constrained QPs (Section VI-B), using a Core i5 -5300U -5th Generation @ 2,3 GHz processor. For each set of executions, the median is shown with a dot. We report base statistics in Table II. The failure rates observed for MOSEK and GUROBI solvers come from the fact, again, they are not able reaching the desired precision ϵ abs = 10 -9 , their outputted solutions being not accurate enough.

C. Non-strongly convex problems

As before, in the spirit of the OSQP [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF], we generate random QPs for which the Hessian H is not strictly positive definite. We can see in Figure 3 that when matrices have 15% of sparsity, OSQP and ProxQP have a similar speed for d ≤ 200. For higher dimension, we observe that ProxQP is approximately 1.8 times faster than OSQP. When sparsity is about 50%, one can see in Figure 4 ProxQP is about 1.8 to 2 times faster for d ≈ 50. When d ≈ 1000, ProxQP is about four times faster.

One can see on Figure 5 that when executed on a more modern computer (see Benchmark setup), performance gains are higher, i.e., for d ≥ 200 our solver is 2 to 3 times faster than the second best solver (it was about about 1.8 times faster before).

We report base statistics in Table III. The failure rates observed for MOSEK and GUROBI solvers come from the fact, again, their solutions not being precise enough.

D. Maros-Mészàros problems

The Maros-Mészàros test set [START_REF] Maros | A repository of convex quadratic programming problems[END_REF] is composed of 138 "hard" QPs. Most of them are sparse and ill-conditioned problems, and they contain up to 90597 variables and 180895 constraints. About 83% of the problems have a sparsity level 4. Solving times for pure inequality constrained QPs with non-strictly positive definite Hessian and 50% of sparsity (Section VI-C), using a Core i5 -5300U -5th Generation @ 2,3 GHz processor. For each set of executions, the median is shown with a dot. lower than 10% for the Hessians H, as well as for the constraint matrices A and C.

As we have currently implemented only a dense version of our solver, we restrict the benchmark to problems whose dimensions (constraints and variables) are below or equal to 10 3 ; that is, a subset of 62 problems of the Maros-Mészàros test set (about 45% of the set), for which two thirds have a sparsity level no larger than 20%.

On this set, we measure the ability of the different solvers to tackle those ill-conditioned problems to high accuracy within the predefined runtime limit (1000 seconds). We report a few statistics for each solver under consideration below, including the failure rate (FR), the sifted geometric means (SGM) as well as the performance profile (see details in, e.g., [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF][START_REF] Hermans | QPALM: A Proximal Augmented Lagrangian Method for Nonconvex Quadratic Programs[END_REF]). a) Failure rates: Failure rates mostly come from saturation effects (time limit reached mostly for OSQP, or outputted solutions being not precise enough for others, in some cases primal or dual infeasibility is detected [39, Section 5.1]). These saturations come at different precision levels. Concerning ProxQP, all proposed solutions are outputted within the time limit of 1000 seconds. However, for three of them, unscaling the scaled proposed solutions through reversed equilibration procedure [39, Section 5.1]) makes them finally close to ϵ abs = 10 -9 but just above the threshold (around ≈ 1.1×10 -9 for example, whereas the scaled solutions satisfy the scaled stopping criterion strictly below 10 -9 ). For other solvers, the situation can be close in some cases, but it has also been observed that saturation effects appears often sooner at higher precision levels (≈ 10 -7 for MOSEK or GUROBI solvers for example).

b) Shifted geometric means: Let t s,p denote the time required for solver s to solve problem p. The shifted geometric mean ts of the runtimes for solver s on problem set P is

ts : = |P| p∈P (t s,p + ζ) -ζ = e 1 |P| p∈P log(ts,p+ζ) -ζ. (37)
The second formulation is used in practice to prevent overflow when computing the product. In this paper, runtimes are expressed in seconds, and a shift of ζ = 10 is used, as in the OSQP API for benchmarks [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF]. As in [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF], we employ the convention that when a solver s fails to solve a problem p (within the time limit of 1000 seconds), the corresponding t s,p is set to the time limit for the computation. We report these experimental results in Table IV for the Maros-Mészàros set. 6. Performance profiles on small to medium-sized Maros-Mészàros problems (see Section VI-D), using a Core i5 -5300U -5th Generation @ 2,3 GHz processor. Target accuracy ϵ abs = 10 -9 (see criterion in (2)) and time limit set to 1000 seconds. The higher the better. c) Performance profiles: We report performance profiles on Figure 6. Performance profiles correspond to the fraction of problems solved as a function of certain runtime (measured in terms of a multiple of the runtime of the fastest solver for that problem). More precisely, if S is the set of solvers tested: r s,p := t s,p min s∈S t s,p ,

denotes the performance ratio for solver s with respect to problem p. The fraction of problems q s (τ ) solved by s within a multiple τ of the best runtime, is then given by q s (τ ) = 1 |P| p∈P,rs,p≤τ 1.

Figure 6 shows the plot of q s (τ ) curves for all the solvers s of our test-bed, when executed in our subset of Maros-Mészàros problems defined above. It turns out that our solver always depicts the best performance profile, with only 3 unsolved problems (see more details on it in paragraph VI-D0a), while other solvers appear to be less efficient and output more unsolved problems: either they reach the 10 3 [s] time-out or are unable to satisfy the desired accuracy.

E. Inverse kinematics and dynamics

Finally, we have benchmarked ProxQP against OSQP (the second best solver of our previous tests bed) on the typical inverse kinematics and dynamics settings. We consider the task of controlling the center of mass of the TALOS humanoid robot [START_REF] Stasse | A new humanoid research platform targeted for industrial applications[END_REF] (which should remain within the convex supporting polygon), keeping the two feet on the ground while also moving both arms to reach two random target placements. We use the Pinocchio library [START_REF] Carpentier | The pinocchio c++ library: A fast and flexible implementation of rigid body dynamics algorithms and their analytical derivatives[END_REF] to compute the kinematic and dynamic quantities associated to each task. ProxQP solves the QP problems in 24 ± 7 µs, while OSQP takes 167 ± 93 µs. For the inverse dynamic task, ProxQP solves it in 25 ± 6 µs while OSQP takes 441 ± 193 µs. Such performances may lead to inverse kinematics or inverse dynamics controllers to be run at very high frequency, an important feature for torquecontrolled robots [START_REF] Herzog | Momentum control with hierarchical inverse dynamics on a torquecontrolled humanoid[END_REF] for instance.

VII. CONCLUSION

In this work, we have proposed a new algorithm for solving generic QPs together with its numerical C++ implementation, motivated by robotic applications. We notably propose to combine the bounded constraint Lagrangian (BCL) globalization strategy [START_REF] Conn | A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds[END_REF] with the solving of (primal-dual) Karush-Kuhn-Tucker conditions associated to QP. The intermediary proximal subproblems are solved via a primal semi-smooth Newton method, potentially initiated at an educated guess that allows to ultimately not requiring any Newton iteration when the iterates get close enough to an optimal solution. We highlight the numerical efficiency of our method on various sets of standard QPs ranging from randomly generated problems strongly inspired by the benchmark suite [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF] to a subset of the Maros-Mészàros problems [START_REF] Maros | A repository of convex quadratic programming problems[END_REF]. Our solver turns out to perform better than modern solvers of the literature on small to medium-size problems (d ≤ 10 3 ), while remaining competitive for larger dimensions. As future work, we plan to extend our solver to also support the case of large-dimensional sparse problems for enabling its use in larger-scale problem setups (d ≫ 10 3 ). This will make ProxQP a suitable, efficient and reliable QP solver to operate on modern optimization problems in robotics and beyond.

We also plan to generalize the exploitation of the primaldual augmented Lagrangian techniques to the context of nonlinear problems (e.g. nonlinear optimal control of robotic systems for solving locomotion and manipulation tasks, etc.) following the approach proposed in [START_REF] Kazdadi | Equality Constrained Differential Dynamic Programming[END_REF] and recently extended in [START_REF] Jallet | Implicit differential dynamic programming[END_REF] and [START_REF] Jallet | Constrained differential dynamic programming: A primal-dual augmented lagrangian approach[END_REF]. From a software perspective, we plan to extend our contribution to account for sparse and matrix-free methods for large QP problems (d > 1000).

  Fig. 2.Solving times for sparse degenerate pure inequality constrained QPs (Section VI-B), using a Core i5 -5300U -5th Generation @ 2,3 GHz processor. For each set of executions, the median is shown with a dot.

Fig. 3 .

 3 Fig.3. Solving times for sparse pure inequality constrained QPs with nonstrictly positive definite Hessian (Section VI-C), using a Core i5 -5300U -5th Generation @ 2,3 GHz processor. For each set of executions, the median is shown with a dot.

  Fig.[START_REF] Conn | A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds[END_REF]. Solving times for pure inequality constrained QPs with non-strictly positive definite Hessian and 50% of sparsity (Section VI-C), using a Core i5 -5300U -5th Generation @ 2,3 GHz processor. For each set of executions, the median is shown with a dot.

Fig. 5 .

 5 Fig.5. Solving times for sparse pure inequality constrained QPs with nonstrictly positive definite Hessian (Section VI-C), using an Intel(R) Core(TM) i7-4790 CPU @ 3.60GHz. For each set of executions, the median is shown with a dot.

  Fig.6. Performance profiles on small to medium-sized Maros-Mészàros problems (see Section VI-D), using a Core i5 -5300U -5th Generation @ 2,3 GHz processor. Target accuracy ϵ abs = 10 -9 (see criterion in (2)) and time limit set to 1000 seconds. The higher the better.

  that for problems of dimensions d ≈ 50 our solver is about 3 times faster than OSQP (the second best solver in this set of experiments).

		osqp proxqp gurobi mosek	qpoases quadprog
	SGM [s] 7.1	1.0	461.2	16065.0 118.5	3.7
	FR (%)	0	0	14	76	0	0
				TABLE II
	SHIFTED GEOMETRIC MEANS (SGM) AND FAILURE RATES (FR) FOR
	SPARSE DEGENERATE PURE INEQUALITY CONSTRAINED QPS
				(SECTION VI-B).

TABLE IV SHIFTED

 IV GEOMETRIC MEANS (SGM) AND FAILURE RATE (FR) FOR MAROS-M ÉSZ ÀROS PROBLEMS (SECTION VI-D).

ACKNOWLEDGMENTS

This work was partly funded by the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute), the European Research Council (grant SEQUOIA 724063) and the Louis Vuitton ENS Chair on Artificial Intelligence.