Stéphane Geller

Christophe Hauser

Frédéric Tronel

Valérie Viet

Triem Tong

Information Flow Control for Intrusion Detection derived from MAC Policy

Keywords: Intrusion detection, Information flow control, Mandatory Access Control

Most of today's MAC implementations can be turned into permissive mode, where no enforcement is performed but alerts are raised instead. This behavior is very close to an anomaly IDS except that the system is configured through a MAC policy. MAC implementations such as SELinux and AppArmor come with a default policy including real life and practical rules ready to be used as is or as a basis for a custom policy. In this paper, we first propose an extension of an IDS based on information flow control. We address issues concerning programs execution and improve its expressiveness in terms of security policy. This extended model can be configured to reach a wide variety of different security goals. Particularly, it allows for information flow checking based on users and/or programs dependent policy rules. Furthermore, suspicious modification of binary programs can be detected to avoid malware execution. We also propose an algorithm for deriving an AppArmor MAC policy into an information flow policy, and thus get the advantage of having a ready to use policy offering good security. We finally show a practical example of deriving such a policy in order to configure our IDS.

I. INTRODUCTION

Over the past years, access control mechanisms in most operating systems have been improved. While traditional discretionary access control (DAC) remains widely used, previous research works on mandatory access control (MAC) have led to implementations in common operating systems, such as Linux, FreeBSD, MacOSX and Windows. Examples include SELinux [START_REF] Vance | Implementing selinux as a linux security module[END_REF], AppArmor [START_REF]Apparmor application security for linux[END_REF], Smack [START_REF] Schaufler | the simplified mandatory access control kernel[END_REF], Tomoyo [START_REF] Harada | Access policy generation system based on process execution history[END_REF]. By using those mechanisms, one can finely control the operations each subject is allowed to perform on the objects of the system. A significant amount of work has been spent on defining default security policies for SELinux and Apparmor, offering rules for a lot of applications. This makes those tools valuable for system administrators, reducing the work needed to set up complex security policies in real life systems.

In the same spirit, several models of information flow control have been proposed to address one persistent weakness of access control models, namely the possibility for users or programs to indirectly and illegally access to pieces of information by collaborating with users who have legal access to it. In this article, we focus on intrusion detection and propose a model of information flow control based on mandatory access control. We perform information flow tracking and consider that the detection of an illegal information flow is an intrusion symptom. As we do not enforce the security policy, an alert is raised in case of intrusion and illegal flows are not forbidden. This behavior is known as permissive mode and is available with most MAC implementations. Our model is inspired from a previous model presented in [START_REF] Viet | Specifying and enforcing a fine-grained information flow policy: Model and experiments[END_REF], with three new major contributions. First, our extended model supervises programs execution. Then, policy expressiveness has been improved by using a generic tag system. It is now possible to specify an information flow policy based either on user rights, programs rights, or both at the same time. Finally, in order to determine a practical security policy for this model, we define an algorithm to derive an information flow policy from an AppArmor MAC policy. The paper is organized as follows. First, we briefly present previous works in the literature, related to MAC and information flow tracking. A formal model for representing information flows is then introduced, and our extended model is presented. Finally, we show how it is possible to derive information flow control from a MAC policy.

II. BACKGROUND AND RELATED WORK

In the following, files, sockets and other resources are referred to as objects while processes are referred to as subjects. Any subject or object containing information is a container of information. Discretionary access control (DAC) is the most commonly used access control mechanism and is the default on UNIX based systems. Access is restricted given the identity and the group of the subject who tries to access to an object. It is said to be discretionary because subjects are allowed to transfer certain permissions to each other at their own discretion. Each subject and object has a set of security attributes, and any operation requires to test that it is conform to the policy. Mandatory access control, at the opposite of discretionary access control, is based on authorization rules (policy) enforced by the operating system. The policy is centrally controlled by a security policy administrator, and users cannot modify it. Regular users cannot declassify information, and it is then possible to verify the policy consistency against a given set of security goals [START_REF] Chen | Information flow query and verification for security policy of security-enhanced linux[END_REF], [START_REF] Herzog Joshua | Information flow in operating systems : Eager formal methods[END_REF]. Therefore, a number of security mechanisms are based on MAC, and MAC is central to our approach. This aspect will be further detailed later in section V-B.

Advances in common operating systems such as Linux and FreeBSD include the introduction of generic access control frameworks, including LSM [14] (Linux Security Modules) and TrustedBSD [START_REF] Watson | The trustedbsd mac framework[END_REF]. LSM has led to several implementations, among which SELinux [START_REF] Vance | Implementing selinux as a linux security module[END_REF], Tomoyo [START_REF] Harada | Access policy generation system based on process execution history[END_REF], Smack [START_REF] Schaufler | the simplified mandatory access control kernel[END_REF] and AppArmor [START_REF]Apparmor application security for linux[END_REF] are the most commonly used. When used in enforcement mode, they block illegal accesses to resources before those can be conducted. When used in permissive mode, their behavior is comparable to a model-based IDS.

SELinux [START_REF] Vance | Implementing selinux as a linux security module[END_REF] is the first security module available in Linux, and it has been designed to implement a flexible MAC mechanism called type enforcement (TE). With type enforcement, all subjects and objects have a type identifier. When accessing an object, a subject must have an authorized type of operation (read, write, etc. . . ) with respect to the object's type, and regardless of its user identity. AppArmor [START_REF]Apparmor application security for linux[END_REF] is a simple MAC implementations available in the Linux kernel as an alternative to SELinux. AppArmor aims at being easier to use and configure than SELinux. It is used by default by Novell in their products and comes with a predefined policy, and a set of generic definitions to ease the creation of new policies. AppArmor will be further detailed later in this paper in section V-B.

Contrary to DAC or MAC systems which ensure security in controlling access to containers of information, information flow control ensures security in preventing illegal information flows. Models of information flow control have been introduced in the eighties by Denning, Biba, Bell and Lapadula [START_REF] Denning | A lattice model of secure information flow[END_REF], [START_REF] Biba | Integrity considerations for secure computer systems[END_REF], [START_REF] Lapadula | Secure computer systems: A mathematical model[END_REF]. These models are the origin of the Multilevel Security (MLS) model. In this model, subjects and objects are labeled with a security level, which represent their sensitivity or clearance. Any information flow from lower-level containers to higher-level containers is illegal. Implementations of MLS models try to precisely observe data manipulations in order to prevent illegal information flows. Flume [START_REF] Krohn | Information flow control for standard os abstractions[END_REF], and Histar [15] are modern implementations of information flow control. Flume is an implementation of distributed information flow control (DIFC) for Linux, acting at the OS level, and using standard OS abstractions (processes, pipes, . . . ). In Flume, processes are confined according to a flow control policy. Histar is an operating system aiming at minimizing the amount of code that must be trusted. It provides a secure operating system using mostly untrusted user-level libraries (the only fully trusted code being the kernel). It uses Asbestos labels on six OS level object types (threads, address spaces, segments, gates, containers and devices).

Blare [17], [START_REF] Viet | Specifying and enforcing a fine-grained information flow policy: Model and experiments[END_REF] is a policy-based intrusion detection system aiming at providing fine grained information flow tracking. Content and containers are distinguished, and information flows are observed using tainting techniques. Contents are the data and containers are physical or logical data storage. Tags are associated to containers, independently describing the content as well as the policy for each container in the system. The information flow policy can either be automatically constructed from a DAC policy or adjusted by an administrator. This model has been implemented in Blare 1 and has proved to be helpful to detect attacks (see [START_REF] Hiet | Policy-based intrusion detection in web applications by monitoring java information flows[END_REF]). Nonetheless the model proposed in [START_REF] Viet | Specifying and enforcing a fine-grained information flow policy: Model and experiments[END_REF] is not aware of execution of programs and processes behavior, and does not 1 Blare is freely available at http://www.rennes.supelec.fr/blare/ take different users into account, which is necessary to ensure a fine observation of information flows. Consequently, in this model and its implementation, the authors were not able to derive a Blare policy from a MAC policy, and illegal flows between processes were ignored. In this article we propose to address this problem and we present in section III an extension of the model introduced in [START_REF] Viet | Specifying and enforcing a fine-grained information flow policy: Model and experiments[END_REF]. We also explain how we can now derive an information flow policy from a MAC policy. As an example, we give a general algorithm to derive a Blare policy from an Apparmor policy.

III. EXTENDED BLARE MODEL

In our extended model, we introduce three kinds of containers : "on-disk" containers such as files are called persistent containers (i.e. long-term storage), "in memory" containers are called volatile containers (memory pages, shared memory, IPC. . . ). Processes are considered as a third kind of containers (even though these are volatile) as they correspond to active subjects as opposed to passive memory. We note C the set of all containers, PC the set of all persistent containers, VC the set of all volatile containers and P the set of all processes. Hence, C = PC ∪ VC ∪ P. As multiple processes can in fact execute the same program (e.g. multiple forks of the apache daemon), we also define Π, the set of all classes of processes. Two processes running the same program are in the same class. In the same manner, we distinguish code of running programs from other passive data. We attach meta-information to each element of information in the system that we want to supervise. We note I the set of all meta-information attached to data (e.g. all the personal data of a user, or the code of the apache daemon stored "on disk" in /usr/bin/apache), and X the set of all meta-information attached to running code (as the "in memory" code of the apache process). In practice, elements of I and elements of X are integers. Finally we also model users and we note U the set of all users.

The information flow policy is divided into three independent parts. The first part (P PC ) defines the authorized combinations of atomic information for the set of persistent containers that we want to supervise, the second part (P U ) defines the allowed combinations of atomic information for the users of the system, and the third part (P Π ) defines the authorized combinations of atomic information for the set of all the classes of processes (each class being attached to the code of a program).

Definition 1 (Information flow policy): An information flow policy is a triplet P = (P PC , P U , P Π ) where2 P PC ⊆ PC × ℘(I ∪ X ), P U ⊆ ℘(I ∪ X ) × U and P Π ⊆ Π × ℘(I ∪ X ).

• A pair (c, a) ∈ P C expresses that the container c is allowed to contain any subset of a • A pair (u, a) ∈ P U expresses that any subset of a can be read or executed by the user u. • A pair (π, a) ∈ P Π expresses that any subset of a can be read or executed by a class of processes π running the same code. Then, we introduce the three following notations P PC (c), P U (u) and P Π (π) whose respective values are {a ∈ ℘

(I ∪ X )|(c, a) ∈ P PC }, {a ∈ ℘(I ∪ X )|(u, a) ∈ P U }, {a ∈ ℘(I ∪ X )|(π, a) ∈ P Π }.
Thus, the definition of the information flow policy is defined for persistent containers, for users, and for classes of processes through sets of rules accurately stating which combinations of atomic information those can receive, and which information are authorized to mix together.

These rules are stored in a distributed fashion : we attach three tags to each container : an information tag, a policy tag and an execute policy tag. Dynamic detection of illegal information flows can then be performed at the container level.

The information tag of a container c ∈ C is a set of elements of I ∪ X describing the content of c (i.e. which elements of information it contains). It is updated everytime an information flow modifies the content of c. Note that because we cannot monitor all the information flows at the OS level without hardware memory tagging [16], the information tag contains an over-approximation of the actual content.

The policy tag of a container c ∈ C is a set of subsets of I ∪ X . It defines all the possible combinations of information that are allowed to flow towards the container. This tag is updated when a modification of the policy affects the container c (in practice, we do not change the policy at runtime).

The execute policy tag of a container c ∈ C is a set of subsets of I ∪ X . It defines all the possible combinations of information that processes running the content of c as code are allowed to read or execute. Thus, it is attached to the code of programs, and as the information tag, it is updated everytime an information flow modifies the content of c (i.e. changing the code of a program changes the policy for running it). Note that it is independent of users rights, those are taken into account at execution time (see IV-2).

IV. INFORMATION FLOW TRACKING IN A LIVE SYSTEM

We denote 0, 1, 2, ..., n, ... the states of the system and we note t i the transition between the states i and i + 1 . These transitions correspond to operations requiring an update of the tags content (fork, execution, creation of objects and information flows).

In the following, we will use the functions itag, ptag and xptag that respectively associate an information tag, a policy tag and an execute policy tag to a container. Definition 2 (Legality of an information flow): An information flow towards a container c happening during the transition t i is legal iff ∃p ∈ ptag(c) i+1 such that itag(c) i+1 ⊆ p. We note itag(c) i+1 ptag(c) i+1 when this condition is verified.

• itag : C → ℘(I ∪ X ), itag i (c)
In the following subsections, we are going to detail accurately the propagation of the tags for each operation.

1) f ork: When a process p forks, a clone q is created. We also clone all of its tags.

2) execution of a program: recall that we consider passive data and running code differently. A data is considered as running code once it is executed by a process through an exec() system call. We note Run : I → X , where Run(d) characterizes the running code out of the execution of a data d at transition t i . After an exec() call, the three tags of the calling process p running the object o on behalf of u are initialized as follows. Its information tag becomes :

itag(p) i+1 := k∈itag(o)i\X {Run(k)}
Note that elements of X from o are discarded, because these meta-information are only used to compute the legality of write and append operations (see IV-6).

Its execute policy tag is initialized to the execute policy tag of o.

xptag(p) i+1 := xptag i (o)

Its policy tag is computed from its execute policy tag and from the legal combinations of information for the user running it :

ptag(p) i+1 = xptag i (o) ⊓ P U (u)}
3) Persistent object creation: When a process p creates a new persistent object o on behalf of u, the new object receives an empty information tag.

itag(o) i+1 = ∅

We associate a policy tag to the new object as follows. The authorized flows towards the created object are the flows composed of atomic informations that are legal for u (i.e. the policy states that u is allowed to access this information).

ptag(o) i+1 = P U (u)}

An execute policy tag is also attached to the new object, and it is set to ⊤ by default : there is no restrictions on code execution as it contains no (executable) information.

xptag(o) i+1 = ⊤ 4) Volatile object creation: When a new volatile object is created, it is assigned an empty information tag as it contains no information yet. It is then updated appropriately when further information flows occur (as it will be detailed in the following).

itag(o) i+1 = ∅

In the same manner, the execute policy tag is initialized to ⊤.

xptag(o) i+1 = ⊤
In this model, we consider that anything is allowed to flow into volatile containers. The legality of information flows involving volatile containers depend on the processes operating on it, i.e. processes can only access to information that matches their policy tag. Thus, the policy tag of volatile containers is initialized to ⊤ so that anything can flow into it.

ptag(o) i+1 = ⊤

Using these tags, we are able to perform dynamic detection of illegal information flows by checking if itag(c) ptag(c) stands everytime a flow occurs.

We classify information flows as read like, write like and append like. When a process p reads from an object o, there is a read like information flow between the containers o and p. In the same way, when a process p replaces the content of an object o, there is a write like information flow. Finally, if the process writes without erasing the existing content of the object, it is an append like information flow. When an information flow occurs, we make an over-approximation of the actual flow, considering all the data that might have flown (i.e. all the information tag from the source). It behaves differently whether a read like, write like or append like information flow happends, and tags are updated accordingly for the process and the object.

5) read like operations: When a read like flow occurs on an object o by a process p, we update the information tag and the execute policy tag of p. Note that, as in the case of the execution of a program, elements of X from the information tag of o are discarded (process execution history is not saved in our model).

itag(p) i+1 := itag(p) i ∪ (itag(o) i \ X ) xptag(p) i+1 := xptag(p) i ⊓ (xptag(o) i )
6) write like operations: If the process p overwrites the content of the object o, we simply replace the information tag and the execute policy tag of the object by those of the process. Note that the elements of X in the information tag are used to check the legality of this write like operation performed by this particular code being executed by p.

itag(o)

i+1 := itag(p) i xptag(o) i+1 := xptag(p) i

7) append like operations:

A process p can also append pieces of information to an object o. In this case, the new information tag of o is the concatenation of both information tags. As for write like operations, the elements of X are used to check the legality of the operation. The execute policy tag is updated so as to match both the execute policy tags of the process and of the object.

itag(o) i+1 := itag(o) i ∪ itag(p) i xptag(o) i+1 := xptag(p) i ⊓ xptag(o) i
The figures 1 and 2 summarize the propagations of the tags in the different cases mentionned above. In previous works, the authors have presented how to automatically derive an information flow policy starting from a DAC policy [START_REF] Hiet | Policy-based intrusion detection in web applications by monitoring java information flows[END_REF] or how administrators can compute an adhoc policy [START_REF] Viet | Specifying and enforcing a fine-grained information flow policy: Model and experiments[END_REF]. In this section, we will present an algorithm to derive a Blare policy from a MAC policy.

A. Initialization of tags starting from an existing policy P = (P PC , P U , P Π ) At initialization time, i.e. the initial state of the system, before we start to track information flows, persistent containers are attached an information tag, a policy tag, and an execute policy tag matching the policy.

1) Initial information tag : a unique meta-information describing the initial content of the container is stored into its information tags. This initial information is considered as being atomic (atomic information are the smallest pieces of information that we are able to distinguish in the system).

2) Initial policy tag: for any persistent container c, the associated policy tag is the set of elements in the policy regarding this container. ∀c ∈ PC, ptag 0 (c) := P PC (c)

3) Initial execute policy tag: for any persistent container c, the associated execute policy tag is the set of elements in the policy regarding the execution of the content of c. We note P class : PC → Π the relation that associates a class of processes to any persistent container c. Any process executing the content of c belongs to this class.

xptag 0 (c) := P Π (P class(c))
If the object does not contain executable code, the corresponding class will not appear in P Π and the xptag will be empty.

B. Deriving an information flow policy from an AppArmor Policy

In the following, an information flow policy (centered on programs) is derived from an AppArmor MAC policy. Such a policy does not specify rules based on users, and thus (P U ) is empty. As Blare monitors information flows, we do not take into consideration access control rules that would not be related to any flow transition. To derive a Blare policy from a set of AppArmor profiles, the following is performed : for each statement in the AppArmor policy, we check whether such a statement is related to a potential information flow, and transform it into a Blare statement if it does. The ability to derive such a Blare policy will be useful for future works in comparing different models in terms of intrusion detection, as each model would be configured with a common security policy.

1) AppArmor: In an AppArmor profile, the permission granted to a program π over a resource o can be one of the following : (r,w,l,m,ix,px,Px,ux,Ux). AppArmor profiles also constrain access to network resources and POSIX capabilities. However, these are access control rules and aren't taken into account in this paper. Instead, possible information flows related to those accesses are captured at another level (i.e. actual illegal flow occurs). Such rules would add false positives and are discarded in our derivation.

Definition 3: An AppArmor policy P is a set of profiles. A profile p ∈ P is a set of rules of the form (o, α) where o is an object and α is a permission. All these rules confine a given program π ∈ Π. Such a profile is defined as : (π, {(o 1 , α 1 ), . . . , (o n , α n )})

For each AppArmor policy statement, if it allows a potential flow between a subject and an object, such as defined in section, we update the Blare tag system accordingly.

2) Algorithm: The following algorithm transforms an Ap-pArmor policy (a set of profiles) into an expression of a Blare policy (set of policy labels on containers). Let P be the set of all the AppArmor profiles in the policy. For any profile p ∈ P , p.container is the container associated to the binary program constrained by p, p.canread() is the list of files on which a read like access is authorized, p.canexec() is the list of executable allowed to be executed, and p.canwrite() is the list of paths where it is allowed to write. T OP represents the set of all atomic information tags in the system (it corresponds to ⊤), inherit(p) : bool returns true if the profile p inherits from its parent's profile and f alse otherwise. unconstrained(p) : bool returns true if the associated program (subject) is unconstrained and f alse if not. Run(I) is defined in section III.

VI. EXAMPLE

The following is an example of intrusion detected by this model, when configured with an information flow policy derived from an AppArmor policy. Here, the security is centered on programs, with no user dependent policy rules. Consider the following AppArmor policy example, where two programs are confined : apache and f tpd. Both own files that the other is not allowed to read. We consider AppArmor being setup in permissive mode, and we compare its behavior to our IDS in terms of detection potential. Now, the following execution sequence takes place (see figure 4). The apache process first reads its configuration file /etc/apache2.conf . Then it reads and interprets /www/index.php, containing a security flaw. Arbitrary code is injected and executed through apache. It introduces a malware in the binary code of /usr/bin/f tpd.

In this first part of the execution, the process running apache is not expected to write into /usr/bin/f tpd : the policy tag of this container is not allowed to receive information by a process running apache. Furthermore, the information apache previously read (and figuring in its information tag) does not belong to the policy tag of /usr/bin/f tpd. This would trigger an alert with both AppArmor (configured in permissive mode) and Blare.

Then, apache runs the modified f tpd. The process running apache is allowed to execute f tpd in the security policy, hence AppArmor would allow this execution. But here, the information tag of f tpd has been modified when the arbitrary code was written into it, and meta-information have been added to it. Those new meta-information do not figure in the policy tag of the process running apache, thus it is not authorized to run f tpd anymore, and this would trigger a second alert for illegal code execution with Blare. This quite simple example reveals one of the major goals of Blare : the security administrator can specify a fine-grained information flow policy including processes behavior. Many real life viruses would trigger alerts in this model as soon as the code of a process is changed or confidential information is moved.

VII. CONCLUSION

In this paper we have presented a model of intrusion detection based on an information flow policy, dynamically checking that it is respected. The policy specifies which information may be combined together and which information the containers are allowed to contain. This model offers high expressiveness since we are able to assign meta-information to any data in the system and to constrain the behavior of programs when those data are involved. The policy expresses restrictions on access to information regardless of where it is located in the system by using a tag system associating meta-information to information containers. We explain how we maintain tags when information flows occur and how we can check if the policy is respected. A central concept of this model is the execution of programs. This model performs dynamic checking at execution time, and is able to detect executions of illegal code or illegal flows of information.

Today's MAC implementations in the Linux kernel come with extensive default security policies. It is possible to set up a policy for the model we propose out of an existing MAC policy. We show how to derive a Blare information flow policy from an AppArmor MAC policy, and give an example of practical use.

In our future works we will focus on an implementation of this model in the Linux kernel as a LSM module. We also aim to further enrich this model on two main aspects. First, a user owning information will be able to declassify it. Second, we will provide a high-level language to specify an information flow policy for Blare.

  is the information tag attached to the container c at state i • ptag : C → ℘(℘(I ∪ X )), ptag i (c) is the policy tag attached to the container c at state i • xptag : PC → ℘(℘(I ∪ X )), xptag i (c) is the execute policy tag attached to the container c at state i The relation ⊓ is defined on the sets of sets as follows : A⊓B = {a∩b|a ∈ A, b ∈ B}. We also introduce the notations ⊤ and ⊥ which are respectively the set of all sets of tags and the set containing the empty set.An information flow towards a container c is legal if and only if the new content of c (characterized by itags(c)) is authorized into c, i.e. it appears in ptags(c).

Fig. 1 .Fig. 2 .

 12 Fig. 1. Diagram with an execution and a write-like operation

  Using the previously introduced algorithm, we can derive a Blare policy and compute its expression on the tag system (the function Run() is written R() in the following table):

	{/usr/bin/apache,
	{(/etc/apache2.conf, w),
	(/etc/apache2.conf, r),
	(/www/index.php,r), (/usr/bin/ftpd, px)}
	}
	{/usr/bin/ftpd,
	{(/etc/ftpd.conf,w), (/etc/ftpd.conf,r),
	(/home/ftpd/data,w)}
	}
	path	itag	ptag	xptag
	/usr/bin/apache	{i1}	{i1}	{R(i1), R(i2), i3, i6}
	/usr/bin/ftpd	{i2}	{i2}	{R(i2), i4}
	/etc/apache2.conf	{i3}	{R(i1), i3, i6}	⊤
	/etc/ftpd.conf	{i4}	{R(i2), i4}	⊤
	/home/ftpd/data	{i5}	{R(i2), i4, i5}	⊤
	/www/index.php	{i6}	{R(i1), i3, i6}	⊤

℘ (A) denotes all the subset of a set A


Mark Schmidt
email: mark.schmidt@inria.fr

Nicolas Le Roux
email: nicolas@le-roux.name

Francis Bach
email: francis.bach@ens.fr

Convergence Rates of Inexact Proximal-Gradient Methods for Convex Optimization

HAL is

Introduction

In recent years the importance of taking advantage of the structure of convex optimization problems has become a topic of intense research in the machine learning community. This is particularly true of techniques for non-smooth optimization, where taking advantage of the structure of non-smooth terms seems to be crucial to obtaining good performance. Proximal-gradient methods and accelerated proximal-gradient methods [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] are among the most important methods for taking advantage of the structure of many of the nonsmooth optimization problems that arise in practice. In particular, these methods address composite optimization problems of the form minimize

x∈R d f (x) := g(x) + h(x), (1)
where g and h are convex functions but only g is smooth. One of the most well-studied instances of this type of problem is 1 -regularized least squares [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF][START_REF] Chen | Atomic decomposition by basis pursuit[END_REF],

minimize x∈R d 1 2 Ax -b 2 + λ x 1 ,
where we use • to denote the standard 2 -norm.

Proximal-gradient methods are an appealing approach for solving these types of non-smooth optimization problems because of their fast theoretical convergence rates and strong practical performance. While classical subgradient methods only achieve an error level on the objective function of O(1/ √ k) after k iterations, proximal-gradient methods have an error of O(1/k) while accelerated proximal-gradient methods futher reduce this to O(1/k 2 ) [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF]. That is, accelerated proximal-gradient methods for non-smooth convex optimization achieve the same optimal convergence rate that accelerated gradient methods achieve for smooth optimization.

Each iteration of a proximal-gradient method requires the calculation of the proximity operator, prox L (y) = arg min

x∈R d L 2 x -y 2 + h(x), (2)
where L is the Lipschitz constant of the gradient of g. We can efficiently compute an analytic solution to this problem for several notable choices of h, including the case of 1 -regularization and disjoint group 1 -regularization [START_REF] Wright | Sparse reconstruction by separable approximation[END_REF][START_REF] Bach | Convex optimization with sparsity-inducing norms[END_REF]. However, in many scenarios the proximity operator may not have an analytic solution, or it may be very expensive to compute this solution exactly. This includes important problems such as total-variation regularization and its generalizations like the graph-guided fused-LASSO [START_REF] Fadili | Total variation projection with first order schemes[END_REF][START_REF] Chen | Graph-structured multi-task regression and an efficient optimization method for general fused Lasso[END_REF], nuclearnorm regularization and other regularizers on the singular values of matrices [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF][START_REF] Ma | Fixed point and Bregman iterative methods for matrix rank minimization[END_REF], and different formulations of overlapping group 1 -regularization with general groups [START_REF] Jacob | Group Lasso with overlap and graph Lasso[END_REF][START_REF] Jenatton | Proximal methods for sparse hierarchical dictionary learning[END_REF]. Despite the difficulty in computing the exact proximity operator for these regularizers, efficient methods have been developed to compute approximate proximity operators in all of these cases; accelerated projected gradient and Newton-like methods that work with a smooth dual problem have been used to compute approximate proximity operators in the context of total-variation regularization [START_REF] Fadili | Total variation projection with first order schemes[END_REF][START_REF] Barbero | Fast Newton-type methods for total variation regularization[END_REF], Krylov subspace methods and low-rank representations have been used to compute approximate proximity operators in the context of nuclear-norm regularization [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF][START_REF] Ma | Fixed point and Bregman iterative methods for matrix rank minimization[END_REF], and variants of Dykstra's algorithm (and related dual methods) have been used to compute approximate proximity operators in the context of overlapping group 1 -regularization [START_REF] Jenatton | Proximal methods for sparse hierarchical dictionary learning[END_REF][START_REF] Liu | Fast overlapping group Lasso[END_REF][START_REF] Schmidt | Convex structure learning in log-linear models: Beyond pairwise potentials[END_REF].

It is known that proximal-gradient methods that use an approximate proximity operator converge under only weak assumptions [START_REF] Patriksson | A unified framework of descent algorithms for nonlinear programs and variational inequalities[END_REF][START_REF] Combettes | Solving monotone inclusions via compositions of nonexpansive averaged operators[END_REF]; we briefly review this and other related work in the next section. However, despite the many recent works showing impressive empirical performance of (accelerated) proximal-gradient methods that use an approximate proximity operator [START_REF] Fadili | Total variation projection with first order schemes[END_REF][START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF][START_REF] Ma | Fixed point and Bregman iterative methods for matrix rank minimization[END_REF][START_REF] Barbero | Fast Newton-type methods for total variation regularization[END_REF][START_REF] Liu | Fast overlapping group Lasso[END_REF][START_REF] Schmidt | Convex structure learning in log-linear models: Beyond pairwise potentials[END_REF], up until recently there was no theoretical analysis on how the error in the calculation of the proximity operator affects the convergence rate of proximal-gradient methods. In this work, we show in several contexts that, provided the error in the proximity operator calculation is controlled in an appropriate way, inexact proximal-gradient strategies achieve the same convergence rates as the corresponding exact methods. In particular, in Section 4 we first consider convex objectives and analyze the inexact proximal-gradient (Proposition 1) and accelerated proximal-gradient (Proposition 2) methods. We then analyze these two algorithms for strongly convex objectives (Proposition 3 and Proposition 4). Note that, in these analyses, we also consider the possibility that there is an error in the calculation of the gradient of g. We then present an experimental comparison of various inexact proximal-gradient strategies in the context of solving a structured sparsity problem (Section 5).

Related Work

The algorithm we shall focus on in this paper is the proximal-gradient method

x k = prox L y k-1 -(1/L)(g (y k-1 ) + e k ) , (3)
where e k is the error in the calculation of the gradient and the proximity problem ( 2) is solved inexactly so that x k has an error of ε k in terms of the proximal objective function [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF].

In the basic proximal-gradient method we choose y k = x k , while in the accelerated proximalgradient method we choose

y k = x k + β k (x k -x k-1 ),
where the sequence {β k } is chosen to accelerate the convergence rate.

There is a substantial amount of work on methods that use an exact proximity operator but have an error in the gradient calculation, corresponding to the special case where ε k = 0 but e k is non-zero. For example, when the e k are independent, zero-mean, and finitevariance random variables, then proximal-gradient methods achieve the (optimal) error level of O(1/ √ k) [START_REF] Duchi | Efficient online and batch learning using forward backward splitting[END_REF][START_REF] Langford | Sparse online learning via truncated gradient[END_REF]. This is different than the scenario we analyze in this paper since we do not assume unbiased nor independent errors but instead consider a sequence of errors converging to 0. This leads to faster convergence rates and makes our analysis applicable to the case of deterministic and even adversarial errors.

Several authors have recently analyzed the case of a fixed deterministic error in the gradient, and shown that accelerated gradient methods achieve the optimal convergence rate up to some accuracy that depends on the fixed error level [START_REF]Smooth optimization with approximate gradient[END_REF][START_REF] Baes | Estimate sequence methods: extensions and approximations[END_REF][START_REF] Devolder | First-order methods of smooth convex optimization with inexact oracle[END_REF], while the earlier work of [START_REF] Nedic | Convergence rate of incremental subgradient algorithms[END_REF] analyzes the gradient method in the context of a fixed error level. This contrasts with our analysis where by allowing the error to change at every iteration we can achieve convergence to the optimal solution. Also, we can tolerate a large error in early iterations when we are far from the solution, which may lead to substantial computational gains. Other authors have analyzed the convergence rate of the gradient and projected-gradient methods with a decreasing sequence of errors [START_REF] Luo | Error bounds and convergence analysis of feasible descent methods: A general approach[END_REF][START_REF] Friedlander | Hybrid deterministic-stochastic methods for data fitting[END_REF] but this analysis does not consider the important class of accelerated gradient methods. In contrast, the analysis of [START_REF] Devolder | First-order methods of smooth convex optimization with inexact oracle[END_REF] allows a decreasing sequence of errors (though convergence rates in this context are not explicitly mentioned) and considers the accelerated projected-gradient method. However, the authors of this work only consider the case of an exact projection step and they assume the availability of an oracle that yields global lower and upper bounds on the function. This non-intuitive oracle leads to a novel analysis of smoothing methods, but also to slower convergence rates than proximal-gradient methods. The analysis of [START_REF] Baes | Estimate sequence methods: extensions and approximations[END_REF] considers errors in both the gradient and projection operators for accelerated projected-gradient methods but requires that the domain of the function is compact. None of these works consider proximal-gradient methods.

In the context of proximal-point algorithms, there is a substantial literature on using inexact proximity operators with a decreasing sequence of errors, dating back to the seminal work of Rockafellar [START_REF] Rockafellar | Monotone operators and the proximal point algorithm[END_REF]. Accelerated proximal-point methods with a decreasing sequence of errors have also been examined, beginning with [START_REF] Güler | New proximal point algorithms for convex minimization[END_REF]. However, unlike proximal-gradient methods where the proximity operator is only computed with respect to the non-smooth function h, proximal-point methods require the calculation of the proximity operator with respect to the full objective function. In the context of composite optimization problems of the form [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], this requires the calculation of the proximity operator with respect to g + h. Since it ignores the structure of the problem, this proximity operator may be as difficult to compute (even approximately) as the minimizer of the original problem.

Convergence of inexact proximal-gradient methods can be established with only weak assumptions on the method used to approximately solve [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF]. For example, we can establish that inexact proximal-gradient methods converge under some closedness assumptions on the mapping induced by the approximate proximity operator, and the assumption that the algorithm used to compute the inexact proximity operator achieves sufficient descent on problem (2) compared to the previous iteration x k-1 [START_REF] Patriksson | A unified framework of descent algorithms for nonlinear programs and variational inequalities[END_REF]. Convergence of inexact proximalgradient methods can also be established under the assumption that the norms of the errors are summable [START_REF] Combettes | Solving monotone inclusions via compositions of nonexpansive averaged operators[END_REF]. However, these prior works did not consider the rate of convergence of inexact proximal-gradient methods, nor did they consider accelerated proximal-gradient methods. Indeed, the authors of [START_REF] Fadili | Total variation projection with first order schemes[END_REF] chose to use the non-accelerated variant of the proximalgradient algorithm since even convergence of the accelerated proximal-gradient method had not been established under an inexact proximity operator.

While preparing the final version of this work, [START_REF] Villa | Accelerated and inexact forward-backward algorithms[END_REF] independently gave an analysis of the accelerated proximal-gradient method with an inexact proximity operator and a decreasing sequence of errors (assuming an exact gradient). Further, their analysis leads to a weaker dependence on the errors than in our Proposition 2. However, while we only assume that the proximal problem can be solved up to a certain accuracy, they make the much stronger assumption that the inexact proximity operator yields an ε k -subdifferential of h [28, Definition 2.1]. Our analysis can be modified to give an improved dependence on the errors under this stronger assumption. In particular, the terms in √ ε i disappear from the expressions of A k , A k and A k . In the case of Propositions 1 and 2, this leads to the optimal convergence rate with a slower decay of ε i . More details may be found after Lemma 2 in the Appendix. More recently, [START_REF] Jiang | An inexact accelerated proximal gradient method for large scale linearly constrained convex SDP[END_REF] gave an alternative analysis of an accelerated proximal-gradient method with an inexact proximity operator and a decreasing sequence of errors (assuming an exact gradient), but under a non-intuitive assumption on the relationship between the approximate solution of the proximal problem and the ε k -subdifferential of h.

Notation and Assumptions

In this work, we assume that the smooth function g in ( 1) is convex and differentiable, and that its gradient g is Lipschitz-continuous with constant L, meaning that for all x and y in R d we have

g (x) -g (y) L x -y .
This is a standard assumption in differentiable optimization, see [30, §2.1.1]. If g is twicedifferentiable, this corresponds to the assumption that the eigenvalues of its Hessian are bounded above by L. In Propositions 3 and 4 only, we will also assume that g is µ-strongly convex (see [30, §2.1.3]), meaning that for all x and y in R d we have

g(y) g(x) + g (x), y -x + µ 2 ||y -x|| 2 .
However, apart from Propositions 3 and 4, we only assume that this holds with µ = 0, which is equivalent to convexity of g.

In contrast to these assumptions on g, we will only assume that h in ( 1) is a lower semicontinuous proper convex function (see [31, §1.2]), but will not assume that h is differentiable or Lipschitz-continuous. This allows h to be any real-valued convex function, but also allows for the possibility that h is an extended real-valued convex function. For example, h could be the indicator function of a convex set, and in this case the proximity operator becomes the projection operator.

We will use x k to denote the parameter vector at iteration k, and x * to denote a minimizer of f . We assume that such an x * exists, but do not assume that it is unique. We use e k to denote the error in the calculation of the gradient at iteration k, and we use ε k to denote the error in the proximal objective function achieved by x k , meaning that

L 2 x k -y 2 + h(x k ) ε k + min x∈R d L 2 x -y 2 + h(x) , (4)
where y = y k-1 -(1/L)(g (y k-1 ) + e k )). Note that the proximal optimization problem (2) is strongly convex and in practice we are often able to obtain such bounds via a duality gap (e.g., see [START_REF] Jenatton | Proximal methods for sparse hierarchical dictionary learning[END_REF] for the case of overlapping group 1 -regularization).

Convergence Rates of Inexact Proximal-Gradient Methods

In this section we present the analysis of the convergence rates of inexact proximal-gradient methods as a function of the sequences of solution accuracies to the proximal problems {ε k }, and the sequences of magnitudes of the errors in the gradient calculations { e k }. We shall use (H) to denote the set of four assumptions which will be made for each proposition:

• g is convex and has L-Lipschitz-continuous gradient;

• h is a lower semi-continuous proper convex function;

• The function f = g + h attains its minimum at a certain x * ∈ R n ;

• x k is an ε k -optimal solution to the proximal problem (2) in the sense of (4).

We first consider the basic proximal-gradient method in the convex case:

Proposition 1 (Basic proximal-gradient method -Convexity) Assume (H) and that we iterate recursion (3) with y k = x k . Then, for all k 1, we have

f 1 k k i=1 x i -f (x * ) L 2k x 0 -x * + 2A k + 2B k 2 , (5)
with

A k = k i=1 e i L + 2ε i L , B k = k i=1 ε i L .
The proof may be found in the Appendix. Note that while we have stated the proposition in terms of the function value achieved by the average of the iterates, it trivially also holds for the iteration that achieves the lowest function value. This result implies that the wellknown O(1/k) convergence rate for the gradient method without errors still holds when both { e k } and { √ ε k } are summable. A sufficient condition to achieve this is for e k and √ ε k to decrease as O(1/k 1+δ ) for any δ > 0. Note that a faster convergence of these two errors will not improve the convergence rate but will yield a better constant factor.

It is interesting to consider what happens if

{ e k } or { √ ε k } is not summable. For instance, if e k and √ ε k decrease as O(1/k), then A k grows as O(log k) (note that B k is always smaller than A k ) and the convergence of the function values is in O log 2 k k
. Finally, a necessary condition to obtain convergence is that the partial sums A k and B k need to be in o( √ k).

We now turn to the case of an accelerated proximal-gradient method. We focus on a basic variant of the algorithm where β k is set to (k -1)/(k + 2) [32, Eq. ( 19) and ( 27)]:

Proposition 2 (Accelerated proximal-gradient method -Convexity) Assume (H) and that we iterate recursion (3) with

y k = x k + k-1 k+2 (x k -x k-1 ).
Then, for all k 1, we have

f (x k ) -f (x * ) 2L (k + 1) 2 x 0 -x * + 2 A k + 2 B k 2 , (6)
with

A k = k i=1 i e i L + 2ε i L , B k = k i=1 i 2 ε i L .
In this case, we require the series {k e k } and {k √ ε k } to be summable to achieve the optimal O(1/k 2 ) rate, which is an (unsurprisingly) stronger constraint than in the basic case. A sufficient condition is for e k and √ ε k to decrease as O(1/k 2+δ ) for any δ > 0.

Note that, as opposed to Proposition 1 that is stated for the average iterate, this bound is for the last iterate x k .

Again, it is interesting to see what happens when the summability assumption is not met.

First, if e k or √ ε k decreases at a rate of O(1/k 2 ), then k( e k + √ e k ) decreases as O(1/k) and A k grows as O(log k) (note that B k is always smaller than A k ), yielding a convergence rate of O log 2 k k 2 for f (x k ) -f (x * )
. Also, and perhaps more interestingly, if e k or √ ε k decreases at a rate of O(1/k), Eq. ( 6) does not guarantee convergence of the function values. More generally, the form of A k and B k indicates that errors have a greater effect on the accelerated method than on the basic method. Hence, as also discussed in [START_REF] Devolder | First-order methods of smooth convex optimization with inexact oracle[END_REF], unlike in the error-free case, the accelerated method may not necessarily be better than the basic method because it is more sensitive to errors in the computation.

In the case where g is strongly convex it is possible to obtain linear convergence rates that depend on the ratio γ = µ/L, as opposed to the sublinear convergence rates discussed above. In particular, we obtain the following convergence rate on the iterates of the basic proximal-gradient method:

Proposition 3 (Basic proximal-gradient method -Strong convexity) Assume (H), that g is µ-strongly convex, and that we iterate recursion (3) with y k = x k . Then, for all k 1, we have:

x k -x * (1 -γ) k ( x 0 -x * + Āk ) , (7)
with

Āk = k i=1 (1 -γ) -i e i L + 2ε i L .
A consequence of this proposition is that we obtain a linear rate of convergence even in the presence of errors, provided that e k and √ ε k decrease linearly to 0. If they do so at a rate of Q < (1 -γ), then the convergence rate of x k -x * is linear with constant (1 -γ), as in the error-free algorithm. If we have Q > (1 -γ), then the convergence of

x k -x * is linear with constant Q . If we have Q = (1 -γ), then x k -x * converges to 0 as O(k (1 -γ) k ) = o [(1 -γ) + δ ] k for all δ > 0.
Finally, we consider the accelerated proximal-gradient algorithm when g is strongly convex. We focus on a basic variant of the algorithm where

β k is set to (1 - √ γ)/(1 + √ γ) [30, §2.2.1]:
Proposition 4 (Accelerated proximal-gradient method -Strong convexity) Assume (H), that g is µ-strongly convex, and that we iterate recursion (3) with

y k = x k + 1- √ γ 1+ √ γ (x k - x k-1 ). Then, for all k 1, we have f (x k ) -f (x * ) (1 - √ γ) k 2(f (x 0 ) -f (x * )) + A k 2 µ + B k 2 , (8)
with

A k = k i=1 e i + 2Lε i (1 - √ γ) -i/2 , B k = k i=1 ε i (1 - √ γ) -i .
Note that while we have stated the result in terms of function values, we obtain an analogous result on the iterates because by strong convexity of f we have

µ 2 ||x k -x * || 2 ≤ f (x k ) -f (x * ).
This proposition implies that we obtain a linear rate of convergence in the presence of errors provided that ||e k || 2 and ε k decrease linearly to 0. If they do so at a rate Q < (1 -√ γ),

then the constant is (1- √ γ), while if Q > (1- √ γ)
then the constant will be Q . Thus, the accelerated inexact proximal-gradient method will have a faster convergence rate than the exact basic proximal-gradient method provided that Q < (1 -γ). Oddly, in our analysis of the strongly convex case, the accelerated method is less sensitive to errors than the basic method. However, unlike the basic method, the accelerated method requires knowing µ in addition to L. If µ is misspecified, then the convergence rate of the accelerated method may be slower than the basic method.

Experiments

We tested the basic inexact proximal-gradient and accelerated proximal-gradient methods on the CUR-like factorization optimization problem introduced in [START_REF] Mairal | Convex and network flow optimization for structured sparsity[END_REF] to approximate a given matrix W , min

X 1 2 W -W XW 2 F + λ row nr i=1 ||X i || p + λ col nc j=1 ||X j || p .
Under an appropriate choice of p, this optimization problem yields a matrix X with sparse rows and sparse columns, meaning that entire rows and columns of the matrix X are set to exactly zero. In [START_REF] Mairal | Convex and network flow optimization for structured sparsity[END_REF], the authors used an accelerated proximal-gradient method and chose p = ∞ since under this choice the proximity operator can be computed exactly. However, this has the undesirable effect that it also encourages all values in the same row (or column) to have the same magnitude. The more natural choice of p = 2 was not explored since in this case there is no known algorithm to exactly compute the proximity operator.

Our experiments focused on the case of p = 2. In this case, it is possible to very quickly compute an approximate proximity operator using the block coordinate descent (BCD) algorithm presented in [START_REF] Jenatton | Proximal methods for sparse hierarchical dictionary learning[END_REF], which is equivalent to the proximal variant of Dykstra's algorithm introduced by [START_REF] Bauschke | A Dykstra-like algorithm for two monotone operators[END_REF]. In our implementation of the BCD method, we alternate between computing the proximity operator with respect to the rows and to the columns. Since the BCD method allows us to compute a duality gap when solving the proximal problem, we can run the method until the duality gap is below a given error threshold ε k to find an x k+1 satisfying (4).

In our experiments, we used the four data sets examined by [START_REF] Mairal | Convex and network flow optimization for structured sparsity[END_REF] 1 and we choose λ row = .01 and λ col = .01, which yielded approximately 25-40% non-zero entries in X (depending on the data set). Rather than assuming we are given the Lipschitz constant L, on the first iteration we set L to 1 and following [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] we double our estimate anytime g

(x k ) > g(y k-1 ) + g (y k-1 ), x k -y k-1 + (L/2)||x k -y k-1 || 2 .
We tested three different ways to terminate the approximate proximal problem, each parameterized by a parameter α:

• ε k = 1/k α : Running the BCD algorithm until the duality gap is below 1/k α .

• ε k = α: Running the BCD algorithm until the duality gap is below α.

• n = α: Running the BCD algorithm for a fixed number of iterations α.

Note that all three strategies lead to global convergence in the case of the basic proximalgradient method, the first two give a convergence rate up to some fixed optimality tolerance, and in this paper we have shown that the first one (for large enough α) yields a convergence rate for an arbitrary optimality tolerance. Note that the iterates produced by the BCD iterations are sparse, so we expected the algorithms to spend the majority of their time solving the proximity problem. Thus, we used the function value against the number of BCD iterations as a measure of performance. We plot the results after 500 BCD iterations for the four data sets for the proximal-gradient method in Figure 1, and the accelerated proximal-gradient method in Figure 2. In these plots, the first column varies α using the choice ε k = 1/k α , the second column varies α using the choice ε k = α, and the third column varies α using the choice n = α. We also include one of the best methods from the first column in the second and third columns as a reference.

In the context of proximal-gradient methods the choice of ε k = 1/k 3 , which is one choice that achieves the fastest convergence rate according to our analysis, gives the best performance across all four data sets. However, in these plots we also see that reasonable performance can be achieved by any of the three strategies above provided that α is chosen carefully. For example, choosing n = 3 or choosing ε k = 10 -6 both give reasonable performance. However, these are only empirical observations for these data sets and they may be ineffective for other data sets or if we change the number of iterations, while we have given theoretical justification for the choice ε k = 1/k 3 .

Similar trends are observed for the case of accelerated proximal-gradient methods, though the choice of ε k = 1/k 3 (which no longer achieves the fastest convergence rate according to our analysis) no longer dominates the other methods in the accelerated setting. For the SRBCT data set the choice ε k = 1/k 4 , which is a choice that achieves the fastest convergence rate up to a poly-logarithmic factor, yields better performance than ε k = 1/k 3 . Interestingly, the only choice that yields the fastest possible convergence rate (ε k = 1/k 5 ) had reasonable performance but did not give the best performance on any data set. This seems to reflect the trade-off between performing inner BCD iterations to achieve a small duality gap and performing outer gradient iterations to decrease the value of f . Also, the constant terms which were not taken into account in the analysis do play an important role here, due to the relatively small number of outer iterations performed.

Discussion

An alternative to inexact proximal methods for solving structured sparsity problems are smoothing methods [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF] and alternating direction methods [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. However, a major disadvantage of both these approaches is that the iterates are not sparse, so they can not take advantage of the sparsity of the problem when running the algorithm. In contrast, the method proposed in this paper has the appealing property that it tends to generate Objective function against number of proximal iterations for the accelerated proximal-gradient method with different strategies for terminating the approximate proximity calculation. From top to bottom we have the 9 Tumors, Brain Tumor1, Leukemia1, and SRBCT data sets. sparse iterates. Further, the accelerated smoothing method only has a convergence rate of O(1/k), and the performance of alternating direction methods is often sensitive to the exact choice of their penalty parameter. On the other hand, while our analysis suggests using a sequence of errors like O(1/k α ) for α large enough, the practical performance of inexact proximal-gradients methods will be sensitive to the exact choice of this sequence.

Although we have illustrated the use of our results in the context of a structured sparsity problem, inexact proximal-gradient methods are also used in other applications such as total-variation [START_REF] Fadili | Total variation projection with first order schemes[END_REF][START_REF] Chen | Graph-structured multi-task regression and an efficient optimization method for general fused Lasso[END_REF] and nuclear-norm [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF][START_REF] Ma | Fixed point and Bregman iterative methods for matrix rank minimization[END_REF] regularization. This work provides a theoretical justification for using inexact proximal-gradient methods in these and other applications, and suggests some guidelines for practioners that do not want to lose the appealing convergence rates of these methods. Further, although our experiments and much of our discussion focus on errors in the calculation of the proximity operator, our analysis also allows for an error in the calculation of the gradient. This may also be useful in a variety of contexts. For example, errors in the calculation of the gradient arise when fitting undirected graphical models and using an iterative method to approximate the gradient of the log-partition function [START_REF] Wainwright | Tree-reweighted belief propagation algorithms and approximate ML estimation by pseudo-moment matching[END_REF]. Other examples include using a reduced set of training examples within kernel methods [START_REF] Kivinen | Online learning with kernels[END_REF] or subsampling to solve semidefinite programming problems [START_REF]Subsampling algorithms for semidefinite programming[END_REF].

In our analysis, we assume that the smoothness constant L is known, but it would be interesting to extend methods for estimating L in the exact case [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF] to the case of inexact algorithms. In the context of accelerated methods for strongly convex optimization, our analysis also assumes that µ is known, and it would be interesting to explore variants that do not make this assumption. We also note that if the basic proximal-gradient method is given knowledge of µ, then our analysis can be modified to obtain a faster linear convergence rate of (1 -γ)/(1 + γ) instead of (1 -γ) for strongly-convex optimization using a step size of 2/(µ + L), see Theorem 2.1.15 of [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF]. Finally, we note that there has been recent interest in inexact proximal Newton-like methods [START_REF] Schmidt | Projected Newton-type methods in machine learning[END_REF], and it would be interesting to analyze the effect of errors on the convergence rates of these methods.

with {S k } an increasing sequence, S 0 u 2 0 and λ i 0 for all i. Then, for all k 1, then

u k 1 2 k i=1 λ i +   S k + 1 2 k i=1 λ i 2   1/2
Proof We prove the result by induction. It is true for k = 0 (by assumption). We assume it is true for k -1, and we denote by v k-1 = max{u 1 , . . . , u k-1 }. From the recursion, we thus get

(u k -λ k /2) 2 S k + λ 2 k 4 + v k-1 k-1 i=1 λ i leading to u k λ k 2 + S k + λ 2 k 4 + v k-1 k-1 i=1 λ i 1/2
and thus

v k max    v k-1 , λ k 2 + S k + λ 2 k 4 + v k-1 k-1 i=1 λ i 1/2    The two terms in the maximum are equal if v 2 k-1 = S k + v k-1 k i=1 λ i , i.e., for v * k-1 = 1 2 k i=1 λ i + S k + 1 2 k i=1 λ i 2 1/2 . If v k-1 v * k-1 , then v k v * k-1 since the two terms in the max are increasing functions of v k-1 . If v k-1 v * k-1 , then v k-1 λ k 2 + S k + λ 2 k 4 + v k-1 k-1 i=1 λ i 1/2
. Hence, v k v k-1 , and the induction hypotheses ensure that the property is satisfied for k.

The following lemma will allow us to characterize the elements of the ε k -subdifferential of h at x k , ∂ ε k h(x k ). As a reminder, the ε-subdifferential of a convex function a at x is the set of vectors y such that a(t) -a(x) y (t -x) -ε for all t.

Lemma 2 If x i is an ε i -optimal solution to the proximal problem (2) in the sense of (4), then there exists f i such that f i 2ε i L and

L y i-1 -x i - 1 L (g (y i-1 ) + e i ) -f i ∈ ∂ ε i h(x i ) .
Proof We first recall some properties of ε-subdifferentials (see, e.g., [START_REF] Bertsekas | Convex Analysis and Optimization[END_REF]Section 4.3] for more details). By definition, x is an ε-minimizer of a convex function a if and only if a(x) inf y∈R n a(y) + ε. This is equivalent to 0 belonging to the ε-subdifferential ∂ ε a(x).

If a = a 1 + a 2 , where both a 1 and a 2 are convex, we have

∂ ε a(x) ⊂ ∂ ε a 1 (x) + ∂ ε a 2 (x). If a 1 (x) = L 2 x -z 2 , then ∂ ε a 1 (x) = y ∈ R n L 2 x -z - y L 2 ε = y ∈ R n , y = Lx -Lz + Lf L 2 f 2 ε .
If a 2 = h and x is an ε-minimizer of a 1 + a 2 , then 0 belongs to

∂ ε a(x). Since ∂ ε a(x) ⊂ ∂ ε a 1 (x) + ∂ ε a 2 (x),
we have that 0 is the sum of an element of ∂ ε a 1 (x) and of an element of ∂ ε h(x). Hence, there is an f such that

Lz -Lx -Lf ∈ ∂ ε h(x) with f 2ε L . ( 9
)
Using z = y i-1 -(1/L)(g (y i-1 ) + e i ) and x = x i , this implies that there exists f i such that

f i 2ε i
L and

L y i-1 -x i - 1 L (g (y i-1 ) + e i ) -f i ∈ ∂ ε i h(x i ) .
In [28, Definition 2.1], Eq. ( 9) is replaced by Lz -Lx ∈ ∂ ε h(x). Hence, their definition of an approximate solution is equivalent to ours but using f = 0. If we replace f i by 0 in the proof of Proposition 2, we get the O(1/k 2 ) convergence rate using any sequence of errors {ε k } necessary to achieve the O(1/k 2 ) rate in [START_REF] Villa | Accelerated and inexact forward-backward algorithms[END_REF]Th. 4.4]. We can also make the same assumption on f in Proposition 1 to achieve the optimal convergence rate with a decay of √ ε k in O 1/k 0.5+δ instead of O 1/k 1+δ .

Basic proximal-gradient method with errors in the convex case

We now give the proof of Proposition of 1.

Proof Since x k is an ε k -optimal solution to the proximal problem (2) in the sense of (4), we can use Lemma 2 to yield that there exists

f k such that f k 2ε k L and L x k-1 -x k - 1 L (g (x k-1 ) + e k ) -f k ∈ ∂ ε i h(x k ) .
We now bound g(x i ) and h(x i ) as follows:

g(x i ) g(x i-1 ) + g (x i-1 ), x i -x i-1 + L 2 x i -x i-1 2
using L-Lipschitz gradient and the convexity of g,

g(x * ) + g (x i-1 ), x i-1 -x * + g (x i-1 ), x i -x i-1 + L 2 x i -x i-1 2
using convexity of g.

Using the ε i -subgradient, we have

h(x i ) h(x * ) -g (x i-1 ) + e i + L(x i + f i -x i-1 ), x i -x * + ε i .
Adding the two together, we get:

f (x i ) = g(x i ) + h(x i ) = f (x * ) + L 2 x i -x i-1 2 -L x i -x i-1 , x i -x * + ε i -e i + Lf i , x i -x * = f (x * ) + L 2 x i -x i-1 , x i -x i-1 -2x i + 2x * + ε i -e i + Lf i , x i -x * = f (x * ) + L 2 x i -x * -(x i-1 -x * ), (x * -x i ) + (x * -x i-1 ) + ε i -e i + Lf i , x i -x * = f (x * ) - L 2 x i -x * 2 + L 2 x i-1 -x * 2 + ε i -e i + Lf i , x i -x * f (x i ) f (x * ) - L 2 x i -x * 2 + L 2 x i-1 -x * 2 + ε i + ( e i + 2Lε i ) • x i -x *
using Cauchy-Schwartz and f i 2ε i L . Moving f (x * ) on the other side and summing from i = 1 to k, we get:

k i=1 [f (x i )-f (x * )] - L 2 x k -x * 2 + L 2 x 0 -x * 2 + k i=1 ε i + k i=1 ( e i + 2Lε i ) • x i -x * , i.e. k i=1 [f (x i ) -f (x * )] + L 2 x k -x * 2 L 2 x 0 -x * 2 + k i=1 ε i + k i=1 ( e i + 2Lε i ) • x i -x * . (10)
Eq. ( 10) has two purposes. The first one is to bound the values of x i -x * using the recursive definition. Once we have a bound on these quantities, we shall be able to bound the function values using only x 0 -x * and the values of the errors.

Bounding x i -x *

We now need to bound the quantities x i -x * in terms of x 0 -x * , e i and ε i . Dropping the first term in Eq. [START_REF] Ma | Fixed point and Bregman iterative methods for matrix rank minimization[END_REF], which is positive due to the optimality of f (x * ), we have:

x k -x * 2 x 0 -x * 2 + 2 L k i=1 ε i + 2 k i=1 e i L + 2ε i L • x i -x *
We now use Lemma 1 (using

S k = x 0 -x * 2 + 2 L k i=1 ε i and λ i = 2 e i L + 2ε i L ) to get x k -x * k i=1 e i L + 2ε i L +   x 0 -x * 2 + 2 L k i=1 ε i + k i=1 e i L + 2ε i L 2   1/2 . Denoting A k = k i=1 e i L + 2ε i L and B k = k i=1 ε i L , we get x k -x * A k + x 0 -x * 2 + 2B k + A 2 k 1/2 .
Since A i and B i are increasing sequences ( e i and ε i being positive), we have for i k

x i -x * A i + x 0 -x * 2 + 2B i + A 2 i 1/2 A k + x 0 -x * 2 + 2B k + A 2 k 1/2 A k + x 0 -x * + 2B k + A k using the positivity of x 0 -x * 2 , B k and A 2 k .

Bounding the function values

Now that we have a common bound for all x i -x * with i k, we can upper-bound the right-hand side of Eq. ( 10) using only terms depending on x 0 -x * , e i and ε i .

Indeed, discarding L 2 x k -x * 2 which is positive, Eq. (10) becomes k i=1 [f (x i ) -f (x * )] L 2 x 0 -x * 2 + LB k + LA k (A k + x 0 -x * + 2B k + A k ) L 2 x 0 -x * 2 + LB k + 2LA 2 k + LA k x 0 -x * + LA k 2B k L 2 x 0 -x * + 2A k + 2B k 2 .
Since f is convex, we get

f 1 k k i=1 x i -f (x * ) 1 k k i=1 [f (x i ) -f (x * )] L 2k x 0 -x * + 2A k + 2B k 2 .

Accelerated proximal-gradient method with errors in the convex case

We now give the proof of Proposition 2.

Proof Defining

θ k = 2/(k + 1) v k = x k-1 + 1 θ k (x k -x k-1 ) ,
we can rewrite the update for y k as

y k = (1 -θ k+1 )x k + θ k+1 v k , because (1 -θ k+1 )x k + θ k+1 v k = (1 - 2 k + 2 )x k + 2 k + 2 [x k-1 + k + 1 2 (x k -x k-1 )] = x k - 2 k + 2 (x k -x k-1 ) + k + 1 k + 2 (x k -x k-1 ) = x k - k -1 k + 2 (x k -x k-1 ) = y k .
Because g is Lipschitz and g is convex, we get for any z that

g(x k ) g(y k-1 ) + g (y k-1 ), x k -y k-1 + L 2 x k -y k-1 2 g(z) + g (y k-1 ), y k-1 -z + g (y k-1 ), x k -y k-1 + L 2 x k -y k-1 2 . Because -[g (y k-1 ) + e k + L(x k + f k -y k-1 )] ∈ ∂ ε k h(x k ), we have for any z that h(x k ) ε k + h(z) + L(y k-1 -x k ) -g (y k-1 ) -e k + Lf k , x k -z = ε k + h(z) + g (y k-1 ), z -x k + L x k -y k-1 , z -x k + e k + Lf k , z -x k
Adding these bounds together gives:

g(x k ) + h(x k ) = f (x k ) ε k + f (z) + L x k -y k-1 , z -x k + L 2 x k -y k-1 2 + e k + Lf k , z -x k Choosing z = θ k x * + (1 -θ k )x k-1 gives f (x k ) ε k + f (θ k x * + (1 -θ)x k-1 ) + L x k -y k-1 , θ k x * + (1 -θ k )x k-1 -x k + L 2 x k -y k-1 2 + e k + Lf k , θ k x * + (1 -θ k )x k-1 -x k ε k + θ k f (x * ) + (1 -θ k )f (x k-1 ) + L x k -y k-1 , θ k x * + (1 -θ k )x k-1 -x k + L 2 x k -y k-1 2 + e k + Lf k , θ k x * + (1 -θ k )x k-1 -x k (11)
using the convexity of f and the fact that

θ k is in [0, 1]. Since θ k x * + (1 -θ k )x k-1 -x k = θ k (x * -v k )
and

x k -y k-1 = θ k v k + (1 -θ k )x k-1 -y k-1 = θ k v k -θ k v k-1 ,
we have

L x k -y k-1 , θ k x * + (1 -θ k )x k-1 -x k = Lθ 2 k v k -v k-1 , x * -v k = -Lθ 2 k v k -x * 2 + Lθ 2 k v k -x * , v k-1 -x * (12) L 2 x k -y k-1 2 = Lθ 2 k 2 v k -v k-1 2 = Lθ 2 k 2 v k -x * 2 + v k-1 -x * 2 -2 v k -x * , v k-1 -x * (13)
e k + Lf k , θ k x * + (1 -θ k )x k-1 -x k = θ k e k + Lf k , x * -v k .
Summing Eq. ( 12) and ( 13), we get

L x k -y k-1 , θ k x * + (1 -θ k )x k-1 -x k + L 2 x k -y k-1 2 = Lθ 2 k 2 v k-1 -x * 2 -v k -x * 2
Moving all function values in Eq. ( 11) to the left-side, we then get

f (x k ) -θ k f (x * ) -(1 -θ k )f (x k-1 ) Lθ 2 k v k-1 -x * 2 -v k -x * 2 + ε k + θ k e k + Lf k , x * -v k .
Reordering the terms and dividing by θ 2 k gives

1 θ 2 k (f (x k )-f (x * ))+ L 2 v k -x * 2 1 -θ k θ 2 k (f (x k-1 )-f (x * ))+ L 2 v k-1 -x * 2 + ε k θ 2 k + 1 θ k e k + Lf k , x * -v k .
Now we use that for all k greater than or equal to 1,

1 -θ k θ 2 k 1 θ 2 k-1
to apply this recursively and obtain

1 θ 2 k (f (x k ) -f (x * )) + L 2 v k -x * 2 1 -θ 0 θ 2 0 (f (x 0 ) -f (x * )) + L 2 v 0 -x * 2 + k i=1 ε i θ 2 i + k i=1 1 θ i ( e i + 2Lε i ) • x * -v i using f i 2ε i
L . Since v 0 = x 0 and θ 0 = 2, we get

f (x k )-f (x * )+ Lθ 2 k 2 v k -x * 2 Lθ 2 k 2 x 0 -x * 2 +θ 2 k k i=1 ε i θ 2 i +θ 2 k k i=1 1 θ i e i + 2Lε i • x * -v i . (14)
As in the previous proof, we will now use Eq. ( 14) to first bound the values of v i -x * then, using these bounds, bound the function values.

Bounding v i -x *

We now need to bound the quantities v i -x * in terms of x 0 -x * , e i and ε i .

v k -x * 2 x 0 -x * 2 + 2 L k i=1 ε i θ 2 i + k i=1 2 θ i e i L + 2ε i L • x * -v i .
Since θ i = 2/(i + 1), 1 θ i = i+1 2 i since i 1. Thus, we have

v k -x * 2 x 0 -x * 2 + 2 L k i=1 i 2 ε i + k i=1 2i e i L + 2ε i L • x * -v i .
From Lemma 1 (using

S k = x 0 -x * 2 + 2 L k i=1 i 2 ε i and λ i = 2i e i L + 2ε i L )
, and

denoting A k = k i=1 i e i L + 2ε i L and B k = k i=1 i 2 ε i L , we get v k -x * A k + x 0 -x * 2 + 2 B k + A 2 k 1/2
.

Since A i and B i are increasing sequences, we also have for i k:

v i -x * A i + x 0 -x * 2 + 2 B i + A 2 i 1/2 x 0 -x * + 2 A i + B 1/2 i √ 2 x 0 -x * + 2 A k + B 1/2 k √ 2 .

Bounding the function values

Dropping

Lθ 2 k 2 v k -x * 2 in Eq. ( 14) (since it is positive), we thus have

f (x k ) -f (x * ) Lθ 2 k 2 x 0 -x * 2 + 2 B k + 2 A k x 0 -x * + 2 A k + 2 B k Lθ 2 k 2 x 0 -x * 2 + 2 B k + 2 A k x 0 -x * + 4 A 2 k + 2 A k 2 B k Lθ 2 k 2 x 0 -x * + 2 A k + 2 B k 2 and 1 θ 2 k (f (x k ) -f (x * )) L 2 x 0 -x * + 2 A k + 2 B k 2 .

Basic proximal-gradient method with errors in the strongly convex case

Below is the proof of Proposition 3

Proof Again, there exists f i such that f i 2ε i

L and

L x i-1 -x i - 1 L (g (x i-1 ) + e i ) -f i ∈ ∂ ε i h(x i ) .
Since x * is optimal, we have that

x * = prox L x * -1 L g (x *
) . We first separate f k , the error in the proximal, from the rest:

x k -x * 2 = prox L x k-1 - 1 L g (x k-1 ) - 1 L e k + f k -prox L x * - 1 L g (x * ) 2 = prox L x k-1 - 1 L g (x k-1 ) - 1 L e k -prox L x * - 1 L g (x * ) 2 + f k 2 + 2 f k , prox L x k-1 - 1 L g (x k-1 ) - 1 L e k -prox L x * - 1 L g (x * ) prox L x k-1 - 1 L g (x k-1 ) - 1 L e k -prox L x * - 1 L g (x * ) 2 + 2ε k L + 2 2ε k L prox L x k-1 - 1 L g (x k-1 ) - 1 L e k -prox L x * - 1 L g (x * )
using Cauchy-Schwartz and f k

2ε k L x k-1 - 1 L g (x k-1 ) - 1 L e k -x * + 1 L g (x * ) 2 + 2ε k L + 2 2ε k L x k-1 - 1 L g (x k-1 ) - 1 L e k -x * + 1 L g (x * )
using the non-expansiveness of the proximal

x k-1 - 1 L g (x k-1 ) - 1 L e k -x * + 1 L g (x * ) 2 + 2ε k L + 2 2ε k L x k-1 -x * - 1 L (g (x k-1 ) -g (x *
)) + e k L using the triangular inequality.

We continue this computation, but now separating e k , the error in the gradient, from the rest:

x k -x * 2 = x k-1 -x * - 1 L (g (x k-1 ) -g (x * )) 2 + e k 2 L 2 - 2 L e k , x k-1 -x * - 1 L (g (x k-1 ) - 1 L g (x * )) + 2ε k L + 2 2ε k L x k-1 -x * - 1 L (g (x k-1 ) -g (x * )) + e k L x k-1 -x * - 1 L (g (x k-1 ) -g (x * )) 2 + e k 2 L 2 + 2 L e k x k-1 -x * - 1 L (g (x k-1 ) -g (x * )) + 2ε k L + 2 2ε k L x k-1 -x * - 1 L (g (x k-1 ) -g (x * )) + e k L using Cauchy-Schwartz x k-1 -x * - 1 L (g (x k-1 ) -g (x )) 2 + e k 2 L 2 + 2ε k L + 2 L 2ε k L e k + 2 e k L + 2 2ε k L x k-1 -x * - 1 L (g (x k-1 ) -g (x * )) .
We now need to bound

x k-1 -x * -1 L (g (x k-1 ) -g (x *
)) to get the final result. We have:

x k-1 -x * - 1 L (g (x k-1 ) -g (x * )) 2 = x k-1 -x * 2 + 1 L 2 g (x k-1 ) -g (x * ) 2 - 2 L g (x k-1 ) -g (x * ), x k-1 -x * x k-1 -x * 2 + 1 L 2 g (x k-1 ) -g (x * ) 2 - 2 L 1 L + µ g (x k-1 ) -g (x * ) 2 + Lµ L + µ x k-1 -x * 2 using theorem 2.1.12 of [30] = (1 - 2µ L + µ ) x k-1 -x * 2 + 1 L 1 L - 2 L + µ g (x k-1 ) -g (x * ) 2 (1 - 2µ L + µ ) x k-1 -x * 2 + µ 2 L ( 1 L - 2 L + µ ) x k-1 -x * 2
using the negativity of 1 L -2 L+µ and the strong convexity of g

= 1 - µ L 2 x k-1 -x * 2 . Thus x k -x * 2 1 - µ L 2 x k-1 -x * 2 + e k 2 L 2 + 2ε k L + 2 L 2ε k L e k + 2 e k L + 2 2ε k L 1 - µ L x k-1 -x * = 1 - µ L x k-1 -x * + e k L + 2ε k L 2 .
Taking the square root of both sides and applying the bound recursively yields

x k -x * 1 - µ L k x 0 -x * + k i=1 1 - µ L k-i e i L + 2ε i L .
6.4 Accelerated proximal-gradient method with errors in the strongly convex case

We now give the proof of Proposition 4.

Proof We have (following [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF])

x k = y k-1 - 1 L g (y k-1 ) .
We define

α 2 k = (1 -α k )α 2 k-1 + µ L α k v k = x k-1 + 1 α k-1 (x k -x k-1 ) θ k = α k -µ L 1 -µ L y k = x k + θ k (v k -x k ) .
If we choose α 0 = √ γ, then this yields

y k = x k + 1 - √ γ 1 + √ γ (x k -x k-1 ) .
We can bound g(x k ) with

g(x k ) g(y k-1 ) + g (y k-1 ), x k -y k-1 + L 2 x k -y k-1 2 using the convexity of g g(z) + g (y k-1 ), y k-1 -z + g (y k-1 ), x k -y k-1 + L 2 x k -y k-1 2 - µ 2 y k-1 -z 2
using the µ-strong convexity of g.

Using Lemma 2, we have that -

[g (y k-1 ) + e k + L(x k + f k -y k-1 )] ∈ ∂ ε k h(x k ). Hence, we have for any z that h(x k ) ε k + h(z) + L(y k-1 -x k ) -g (y k-1 ) -e k -Lf k , x k -z = ε k + h(z) + g (y k-1 ), z -x k + L x k -y k-1 , z -x k + e k + Lf k , z -x k
Adding these two bounds, we get for any z

f (x k ) ε k + f (z) + L x k -y k-1 , z -x k + L 2 x k -y k-1 2 - µ 2 y k-1 -z 2 + e k + Lf k , z -x k . Using z = α k-1 x * + (1 -α k-1 )x k-1 , we get f (x k ) ε k + f (α k-1 x * + (1 -α k-1 )x k-1 ) + L x k -y k-1 , α k-1 x * + (1 -α k-1 )x k-1 -x k + L 2 x k -y k-1 2 - µ 2 y k-1 -α k-1 x * -(1 -α k-1 )x k-1 2 + e k + Lf k , α k-1 x * + (1 -α k-1 )x k-1 -x k ε k + α k-1 f (x * ) + (1 -α k-1 )f (x k-1 ) + L x k -y k-1 , α k-1 x * + (1 -α k-1 )x k-1 -x k + L 2 x k -y k-1 2 - µ 2 y k-1 -α k-1 x * -(1 -α k-1 )x k-1 2 - µ 2 α k-1 (1 -α k-1 ) x * -x k-1 2 + e k + Lf k , α k-1 x * + (1 -α k-1 )x k-1 -x k
using the µ-strong convexity of f .

We can replace x k -y k-1 using

x k -y k-1 = x k -x k-1 -θ k-1 (v k-1 -x k-1 ) = θ k-1 x k-1 + θ k-1 α k-1 (x k -x k-1 ) + 1 - θ k-1 α k-1 (x k -x k-1 ) -θ k-1 v k-1 = θ k-1 (v k -v k-1 ) + 1 - θ k-1 α k-1 (x k -x k-1
) .

We also have

(1 -α k-1 )x k-1 -x k = -α k-1 v k α k-1 x * + (1 -α k-1 )x k-1 -x k = α k-1 (x * -v k ) , and
y k-1 -α k-1 x * -(1 -α k-1 )x k-1 = y k-1 -α k-1 (x * -v k ) -x k = α k-1 (v k -x * ) -θ k-1 (v k -v k-1 ) -1 - θ k-1 α k-1 (x k -x k-1 )
Thus,

f (x k ) ε k + α k-1 f (x * ) + (1 -α k-1 )f (x k-1 ) -Lθ k-1 α k-1 v k -v k-1 , v k -x * -L(α k-1 -θ k-1 ) x k -x k-1 , v k -x * + Lθ 2 k-1 2 v k -v k-1 2 + L 1 - θ k-1 α k-1 2 2 x k -x k-1 2 + Lθ k-1 1 - θ k-1 α k-1 v k -v k-1 , x k -x k-1 - µα 2 k-1 2 v k -x * 2 - µθ 2 k-1 2 v k -v k-1 2 - µ 1 - θ k-1 α k-1 2 2 x k -x k-1 2 + µα k-1 θ k-1 v k -x * , v k -v k-1 + µ(α k-1 -θ k-1 ) v k -x * , x k -x k-1 -µθ k-1 1 - θ k-1 α k-1 v k -v k-1 , x k -x k-1 - µ 2 α k-1 (1 -α k-1 ) x * -x k-1 2 + α k-1 e k + Lf k , x * -v k .
To avoid unnecessary clutter, we shall denote E k the additional term induced by the errors, i.e.

E k = ε k + α k-1 e k + Lf k , x * -v k .
Before reordering the terms together, we shall also replace all instances of v k -v k-1 with v k -x * -(v k-1 -x * ):

f (x k ) E k + α k-1 f (x * ) + (1 -α k-1 )f (x k-1 ) -Lθ k-1 α k-1 v k -x * 2 + Lθ k-1 α k-1 v k-1 -x * , v k -x * -L(α k-1 -θ k-1 ) x k -x k-1 , v k -x * + (L -µ)θ 2 k-1 2 v k -x * 2 + (L -µ)θ 2 k-1 2 v k-1 -x * 2 -(L -µ)θ 2 k-1 v k -x * , v k-1 -x * + (L -µ) 1 - θ k-1 α k-1 2 2 x k -x k-1 2 + Lθ k-1 1 - θ k-1 α k-1 v k -x * , x k -x k-1 -Lθ k-1 1 - θ k-1 α k-1 v k-1 -x * , x k -x k-1 - µα 2 k-1 2 v k -x * 2 + µα k-1 θ k-1 v k -x * 2 -µα k-1 θ k-1 v k -x * , v k-1 -x * + µ(α k-1 -θ k-1 ) v k -x * , x k -x k-1 -µθ k-1 1 - θ k-1 α k-1 v k -x * , x k -x k-1 + µθ k-1 1 - θ k-1 α k-1 v k-1 -x * , x k -x k-1 - µ 2 α k-1 (1 -α k-1 ) x * -x k-1 2 .
With a bit of well-needed cleaning, this becomes

f (x k ) E k + α k-1 f (x * ) + (1 -α k-1 )f (x k-1 ) + L -µ 2 (θ k-1 -α k-1 ) 2 - Lα 2 k-1 2 v k -x * 2 + (L -µ)θ 2 k-1 2 v k-1 -x * 2 + (L -µ)θ k-1 (α k-1 -θ k-1 ) v k-1 -x * , v k -x * + (L -µ)(θ k-1 -α k-1 ) 1 - θ k-1 α k-1 x k -x k-1 , v k -x * -(L -µ)θ k-1 1 - θ k-1 α k-1 v k-1 -x * , x k -x k-1 + (L -µ) 1 - θ k-1 α k-1 2 2 x k -x k-1 2 - µ 2 α k-1 (1 -α k-1 ) x * -x k-1 2 .
We can rewrite x k -x k-1 using

x k -x k-1 = α k-1 (v k -x k-1 ) = α k-1 (v k -x * ) -α k-1 (x k-1 -x * ) .
We may now compute the coefficients for the following terms: v k -x * 2 , v k-1 -x * 2 , v k-1 -x * , v k -x * , x k-1 -x * , v k -x * , x k-1 -x * , v k-1 -x * and x * -x k-1 2 .

For v k -x * 2 , we have

L -µ 2 (θ k-1 -α k-1 ) 2 - Lα 2 k-1 2 v k -x * 2 term -(L -µ)(θ k-1 -α k-1 ) 2 x k -x k-1 , v k -x * term + (L -µ)(θ k-1 -α k-1 ) 2 2 x k -x k-1 2 term = - Lα 2 k-1 2 .
For v k-1 -x * 2 , there is only one term and we keep

(L-µ)θ 2 k-1 2
.

For v k-1 -x * , v k -x * , we get

(L -µ)θ k-1 (α k-1 -θ k-1 ) v k-1 -x * , v k -x * term -(L -µ)θ k-1 (α k-1 -θ k-1 ) v k-1 -x * , x k -x k-1 term = 0 . For x k-1 -x * , v k -x * , we get (L -µ)(θ k-1 -α k-1 ) 2 x k -x k-1 , v k -x * term -(L -µ)(θ k-1 -α k-1 ) 2
x k -x k-1 2 term = 0 .

Reordering the terms, we have

f (x k )-f (x * )+ Lα 2 k-1 2 v k -x * 2 (1-α k-1 ) (f (x k-1 ) -f (x * ))+ (Lα k-1 -µ)α k-1 2 v k-1 -x * 2 +E k .
(15) We can rewrite Eq. (15) as

f (x k ) -f (x * ) + Lα 2 k-1 2 v k -x * 2 (1 -α k-1 ) f (x k-1 ) -f (x * ) + Lα 2 k-2 2 v k-1 -x * 2 + Lα 2 k-1 -µα k-1 -(1 -α k-1 )Lα 2 k-2 2 v k-1 -x * 2 + E k .
Using

α k = µ L
and denoting

δ k = f (x k ) -f (x * ) + µ 2 v k -x * 2 , (16)
we get the following recursion:

δ k 1 - µ L δ k-1 + E k . (17)
Applying this relationship recursively, we get

δ k 1 - µ L k δ 0 + k t=1 E t 1 - µ L k-t . (18)
Since E k = ε k + α k-1 e k + Lf k , x * -v k , we can bound it by

E k ε k + µ L e k + 2Lε k v k -x * (19)
using f k 2ε k L . Plugging Eq. ( 19) into Eq. ( 18), we get

δ k 1 - µ L k δ 0 + k t=1 ε t + µ L e t + 2Lε t v t -x * 1 - µ L -t . (20)
Again, we shall use Eq. ( 18) to first bound the values of v i -x * , then the function values themselves.

Figure 1 :

 1 Figure1: Objective function against number of proximal iterations for the proximal-gradient method with different strategies for terminating the approximate proximity calculation. From top to bottom we have the 9 Tumors, Brain Tumor1, Leukemia1, and SRBCT data sets.

Figure 2 :

 2 Figure2: Objective function against number of proximal iterations for the accelerated proximal-gradient method with different strategies for terminating the approximate proximity calculation. From top to bottom we have the 9 Tumors, Brain Tumor1, Leukemia1, and SRBCT data sets.

The datasets are freely available at http://www.gems-system.org.

Acknowledgements

Mark Schmidt, Nicolas Le Roux, and Francis Bach are supported by the European Research Council (SIERRA-ERC-239993).

Appendix: Proofs of the propositions

We first prove a lemma which will be used for the propositions.

Lemma 1 Assume that the nonnegative sequence {u k } satisfies the following recursion for all k 1:

.

Hence, we have

the last two lines allowing us to complete the square. We may now factor it to get

Discarding the term depending on x k-1 -x * and regrouping the terms depending on v k-1x * 2 , we have

We will now use Lemma 1 to bound the value of v k -x * . Since v k -x * 2 is bounded by 2δ k µ (using Eq. ( 16)), we can use Eq. ( 18) to get

Multiplying both sides by 1

Using Lemma 1 with

, we get

with

Since B t is an increasing sequence and the λ t are positive, we have for i k

Hence,

Bounding the function values

Denoting

, we obtain after plugging Eq. ( 23) into Eq. ( 20)

Using the L-Lipschitz gradient of f and the fact that v 0 = x 0 , we have

Hence, discarding the term v k -x * 2 of δ k , we have