Learning With Subquadratic Regularization : A Primal-Dual Approach

1 Introduction

Structured sparse regularizers [Kyrillidis, 2016; Bach et al., 2011; Jenatton et al., 2011b; Micchelli et al., 2013] have emerged as efficient and versatile tools to add prior knowledge to estimation problems arising in domains such as computer vision [Mairal et al., 2014], bioinformatics [Obozinski et al., 2011], neural imaging [Jenatton et al., 2011a] among many others. They typically lead to convex optimization problems of the form

where X ∈ Rn×d is the data matrix, w ↦→ F (Xw) is the convex data-fitting term, the norm Ω is the regularizer, and λ > 0 is a balancing hyper-parameter. We refer to F : Rn → R as the loss function. Popular loss functions include the square loss used in least-squares regression [Tibshirani, 1994], the hinge loss used in support vector machines (SVM) [Vapnik, 2000], or the logistic loss used in logistic regression [Lee et al., 2006]. Many algorithms exist if the proximal operator for Ω (denoted as proxΩ) is known and efficient to compute [Bach et al., 2011; Parikh and Boyd, 2014]. Depending on the loss function F, we may use forward-backward splitting methods such as FISTA [Beck and Teboulle, 2009; Nesterov, 2007] when F is smooth; ADMM [Parikh and Boyd, 2014], when one knows the proximal operator of F ◦X : w ↦→ F (Xw); or Chambolle-Pock algorithm [Chambolle and Pock, 2011], when only proxF is known.

∗The author is currently affiliated with LinkedIn, while the work was carried out when the affiliation was Indian Institute of Science.

In this paper, we focus on norms Ω which are subquadratic [Bach et al., 2011], expressed as follows:

where Γ : Rd+ ↦→ R is convex and positively homogeneous. Note that the norm (2) includes as special cases the popular examples such as the ℓ1-norm [Tibshirani, 1994] and grouped ℓ1,p-norms for p ∈ [1, 2] [Jenatton et al., 2011b; Kloft et al., 2011]. For many such norms of the form (2), neither Ω nor proxΩ is easily evaluated, where proxΓ may be computable with an existing efficient algorithm. Examples include various norms derived from tree and graph structured constraints [Micchelli et al., 2013; Baldassarre et al., 2012], box-structured constraints [McDonald et al., 2016] and norms obtained as convex relaxations of combinatorial penalties [Obozinski and Bach, 2016; Sankaran et al., 2017], among many others. Thus, to solve (1) with Ω of the form (2), we may not be able to easily extend the aforementioned techniques which rely on efficient computation of proxΩ. Hence, we use the relation (2) and reformulate (1) into the following optimization problem in variables w and η, which is the main focus of this paper:

To solve (3), when proxF◦X and proxΓ are easy to compute, one can devise an alternating optimization routine (minimizing w.r.t w and η alternatively until convergence). But this has convergence issues in general because the objective function is not smooth around vectors with zero elements, and does not offer a convergence rate. Alternatively, under the same assumptions, we can pose (3) as a smooth problem in η alone, which can be solved using FISTA, but suffers from expensive per-step computation. The only generic (i.e., with no strong assumptions on F) first-order algorithm is subgradient descent, which is very slow with a rate of O(1/√T ) after T iterations. Thus, we are in need of efficient algorithms to solve problem (3). We make the following contributions:

• We investigate primal-dual algorithms for the problem (3). Though it is complicated to derive saddle point problems directly from (3), we propose a conic reformulation of the problem (3) in Section 5.1, which opens the possibilities for applying efficient primal-dual procedures. Also, the cone constraints derived are separable over the set of variables involved, and have closed form solutions for projection, which is used by the primaldual algorithms.
• When proxF◦X and proxΓ are easy to compute, we propose in Section 5.2, a primal-dual algorithm ADMM-η to solve problem (3) which converges at the rate of O(1/T) in T iterations (see Theorem 1). To achieve an ϵ-approximate solution with respect to the duality gap, ADMM-η takes O((d+ cF◦X + cΓ) /ϵ), where cF◦X and cΓ denote the complexity to evaluate the proximal operator for F ◦X and Γ respectively. In general, cF◦X may be high because of the dependence on X. However for popular loss functions such as the squared loss and hinge loss, this computation turns out to be simple (see Table 1).
• When only proxF and proxΓ are easy to compute, we propose a generic first-order primal-dual algorithm CP-η in Section 5.3, which also converges at a rate of O(1/T ) (See Theorem 2). To guarantee an ϵ-approximate solution, CP-η takes O((nd+ cF + cΓ) /ϵ) operations, where cF denotes the complexity to evaluate the proximal operator for F. In situations where cF◦X is high, CP-η serves as a effective alternative to ADMM-η, while also being the most generic and efficient solution for (3).
• In Section 6, we consider an example of subquadratic norm (See Example 1) which encourages tree sparsity. As studied in [Obozinski and Bach, 2016], proxΓ can be evaluated in O(d log d) time, whereas proxΩ requires O(d2 ) computations. We illustrate that ADMM-η and CP-η outperform existing benchmarks on this chosen example.

Notations. 1d (resp. 0d) denotes a vector in Rd with all 1's (resp. 0's). Given z ∈ Rd , we denote by D(z) the diagonal matrix formed with zi in the (i, i)th entry, and define Id = D(1d). Given a set A, define IA(a) = 0 if a ∈ A, and +∞otherwise. Given f : Rn → R, the Fenchel dual of f defined as f ∗ (α) = supx∈Rn x⊤α-f(x). The proximal operator of a function f at z defined as proxf (z) = argminx 1 2 ‖xz‖2 + f(x), and we define proxτf (z) = proxτf(z). We define proxF◦X(z) = argminw 12 ‖w-z‖2 +F(Xw). Following [Nesterov, 2013], we define a function f as L-smooth if ∀x, y ∈ Rn , f(y) ≤ f(x) + ∇f(x)⊤(y - x) + L2 ‖y - x‖22 .

We define f as µ-strongly convex if ∀x, y ∈ Rn , f(y) ≥ f (x) + ∇f (x)⊤(y - x) + µ2 ‖y - x‖22 .

2 Subquadratic Norms

Many existing studies have illustrated benefit of using subquadratic norms for structured sparsity in many applications: compressed sensing and image denoising [Baldassarre et al., 2012], multi task learning [McDonald et al., 2016], banded covariance matrix estimation [Yan and Bien, 2015]. Next, we discuss specific examples where proxΓ is easy to compute.

Example 1. Tree-Structured H-norms [Baldassarre et al., 2012; Micchelli et al., 2013]. Consider a rooted tree with d nodes represented by the edge matrix A ∈ Rd×d , where Aei = 1,Aej = -1 if eth edge links i to j. Define H = {η ∈ Rd+|Aη ∈ Rd+}, and let Γ(η) = η⊤1d + 1H(η). With this choice of Γ, let us denote the norm (2) as ΩH. By enforcing the tree structured prior on η, ΩH thus encourages w to be tree structured. The proximal operator for ΩH may be computed approximately using Picard iterates [Baldassarre et al., 2012], which also lacks convergence rates. Whereas, we can compute proxΓ easily since the projection onto H is computed easily through a pooladjacent-violators (PAV) algorithm [Pardalos and Xue, 1999; Best and Chakravarti, 1990] in O(d log d) time.

Example 2. Convex Relaxation of Combinatorial Penalties [Bach, 2010; Bach, 2011b; Obozinski and Bach, 2016]. Denote V = {1, . . . , d}. Given a tree structured ordering over V , define Gj = {j ∪ Dj}, where Dj denotes the descendants of j. We may arrive at a subquadratic regualarizer which encourages tree structures through the steps discussed next. Consider a set function S(A) = ∑Gj∈G 1Gj∩A̸=∅. S is submodular [Bach, 2011a], and its corresponding Lov´asz extension denoted by Γ(η) equals ∑j ‖ηGj ‖∞, where ηGj is the restriction of η to the coordinates given by Gj. With this choice of Γ, we denote the norm (2) as ΩS2 (w). proxΩ is more expensive in this case: can be computed using a divide-andconquer strategy involving a sequence of submodular function minimization (SFM), whose complexity is d(O(SFM)). Whereas, proxΓ can be computed in time O(dh), where h is the depth of the tree[Jenatton et al., 2011b].

Comparing Examples 1 and 2, note that ΩS2 (w) = ΩH(w) [Obozinski and Bach, 2016], illustrating that Γ need not be unique for a given subquadratic norm Ω. While these existing works on subquadratic norms illustrate the benefits in terms of applicability, they do not focus on computational efficiency, which is the focus of this paper.

3 Existing Algorithms

FISTA-η. Defining a smooth function J(η) as follows:

problem (3) is equivalent to minη≥0 J(η) + λ2 Γ(η) on which we can use (accelerated) proximal gradient descent [Beck and Teboulle, 2009] with J as the smooth component, and

Γ(η) as the non-smooth component, whose proximal operator will be needed. The key is the computation of the gradient of J, for which we solve (4). With the proper change of variable, the minimizer of (4) is equivalent to w = D(η)1/2prox1/λ F◦XD(η)1/2 (0). When the Lipshitz constant of the gradient of J is not known, backtracking is used to find the stepsize [Scheinberg et al., 2014], which increases the number of times the gradient of J is calculated. We denote this algorithm as FISTA-η. Note that the per-iteration cost may overweigh the benefit of the fast convergence rate O(1/T 2 ) of FISTA-η.

Alt-η. A simple algorithm is to alternatively solve for w and η at each iteration. Given η, to solve for w, we get an ℓ2-regularized problem on w, whose solution is given as w = D(η)1/2prox1/λ F◦XD(η)1/2 (0). Given w, solving for η is equivalent to evaluating the norm Ω(w) and returning the minimizer η. We refer to this algorithm as Alt-η. Though lacking in convergence guarantees and having issues with convergence when η has components close to 0, this algorithm is used in practice with good performances. Initializing η far from 0 is recommended or making sure that η stays away from 0 by adding a penalty ε ∑dj=1 η-1 j [Canu and Grandvalet, 1999; Daubechies et al., 2010; Argyriou et al., 2008].

4 Primal-Dual Algorithms

Before discussing our proposed algorithms for solving (3), we briefly discuss the generic primal-dual setup we work on in this paper. Following [Chambolle and Pock, 2011], for finitedimensional vector spaces U and V, we consider a generic primal problem of the form

where K : U → V is a linear operator from U to V with the operator norm ‖K‖ = max{‖Ku‖, u ∈ U, ‖u‖ ≤ 1}, H : V → R ∪ {+∞} and G : U → R ∪ {+∞} are proper convex functions, whose proximal operators are easily computable. This leads to the following saddle-point problem.

Additionally when H is (1/γ)-smooth, H∗ is γ-strongly convex. The primal-dual algorithm (CP) [Chambolle and Pock, 2011] for (6) by is given in Algorithm 1, which guarantees that there exists positive R1,R2 such that after T steps the following bounds hold true:

the second inequality being valid when F ∗ is γ-stronglyconvex. The key features of Algorithm 1 are as follows:

1. It requires as input proxG and proxH∗, initial primal and dual step-sizes τ (0) , σ(0) satisfying τ (0)σ(0) ‖K‖2 ⩽ 1.
2. It accesses K through only matrix-vector products, and hence strictly first-order in nature.

Algorithm 1 CP [Chambolle and Pock, 2011]

Require: K, proxG, proxH∗ , u(0) , v(0) , T

• 1: Choose τ (0) , σ(0) > 0 such that τ (0)σ(0) ‖K‖2 ⩽ 1.
• 2: Initialize ¯u(0) = u(0)
• 3: for k = 0, 1, 2, ..., T -1 do
• 4: v(k+1) = proxσ(k)H∗(v(k) +σK˜u(k) )
• 5: u(k+1) = proxτ (k)G(u(k) - τK⊤v(k+1) )
• 6: θ(k) = 1 1+2γτ (k) , τ (k) = θ(k) τ (k) , σk = τ (k) θ(k)
• 7: ˜u(k+1) = u(k+1) + θ(k) (u(k+1) -u(k) )
• 8: end for
• 9: Return (u(1:T) , v(1:T) ).
• The cost per iteration is cu + cv + cost(proxG) + cost(prox∗H ), where cu and cv are the cost to compute Ku and K⊤v respectively for vectors u ∈ U, v ∈ V.
• coincides with ADMM, when K = I.

4.1 Step-Size Selection

Given (u, v) ∈ U ×V, referring to Theorem 1 of [Chambolle and Pock, 2011], when γ = 0, the following bound holds for all σ(0) = σ, τ (0) = τ satisfying στ ‖K‖2 ≤ 1.

Let (u∗ , v∗ ) be an optimal solution of (6). Using the constraint στ‖K‖2 ≤ 1, we optimize the upper bound in (9) for the choice (u, v) = (u∗ , v∗ ) resulting in the value σ = 1 ‖K‖ ‖v∗-v(0) ‖2 ‖u∗-u(0) ‖2 . Since it is not possible to calculate σ in practice because we do not have access to u∗ or v∗ , we will compute rough estimates using the information we have. In the experiments, we derive the step-sizes for the squared loss F(z) = 1 2n ‖z-y‖22 and the ℓ1 norm Γ(η) = η⊤1d. Note that the obtained step-sizes work well enough in our experiments, but that by an additional tuning, they could be improved.

5 Primal-Dual Formulation for Subquadratic Norms

We first derive a conic reformulation of (3), which gets rid of the ratio term w2 j 2ηj , enabling application of Algorithm 1.

5.1 Conic-Constrained Primal Reformulation

Note that the term w2 j 2ηj may be written as w2 j 2ηj =mintj tj such that tj ⩾ 0,w2 j ⩽ 2tjηj. Now, Eq. (3) is equivalent to:

where C = {(a, b, c) ∈ R3 , a2 ⩽ 2bc, b ⩾ 0, c ⩾ 0} is the rotated second order-cone. Note that the cone constraint is separable across the sets of variables (wj , ηj , tj ) and thus can be handled separately. The cone C is self-dual and the proximal operator for ∑dj=1 IC(wj, ηj, tj), which will be the key to designing the primal-dual algorithms, involves computing the orthogonal projection onto C using a simple closed form expression. Thus, (12) makes it easy for deriving saddle-point problems of the form (6), which can be efficiently solved using Algorithm (1). In the next subsections, we consider reformulations of (12) under cases - (a) proxF◦X is easy to compute, (b) only proxF is computable easily.

5.2 ADMM-η: when proxF◦X is Computable

When proxF◦X is computable easily, we propose a primaldual algorithm ADMM-η to solve (12). Algorithm 2 lists the steps, which accepts proxΓ and proxF◦X inputs. The follow-

Algorithm 2 ADMM-η

Require: X, proxΓ, proxF◦X, k

• 1: Initialize u(0) , v(0) ∈ R3d , K = -I3d. (Ref. (15)).
• 2: Define functions G, H∗ as in (16).
• 3: (u(k) , v(k) ) = CP(K, proxG, proxH∗, u(0) , v(0) , k)
• 4: w(k) = u(k) (1 : d), η(k) = u(k) (d + 1 : 2d).
• 5: Return w(k) , η(k) .

ing result proves that Algorithm 2 results in a rate of O(1/T ) for the problem (3).

Theorem 1. Let F : Rn → R and Γ : Rn+ → R be convex, and let Ψ(w, η) be as defined in (3) with (w∗ , η∗ ) = argminw,η Ψ(w, η), and (u(k) = [w(k) η(k) t(k) ]⊤ , v(k) = [δ(k)β(k)γ (k) ]⊤ )k≥0 be defined as in Algorithm 2. Then, there exists a constant R > 0 such that after T iterations

Proof. (Sketch) First, we can rewrite the cone constraint IC in (12) in terms of its Fenchel dual, which again is a similar cone constraint on the dual variables, arriving at the following saddle-point problem equivalent to (12):

-

We can equate the above problem to (6) through the mapping of the variables u, v and functions G, H∗ as given below.

Since the function G is separable in terms of the primal variables w, η, t, proxG simply reduces to computing the proximal operators independently on each group of variables.

proxG and prox∗H may be computed straightforward using proxF◦X and proxΓ. To evaluate proxG, the cost is usually dominated by the complexity to solve proxF◦X. Now result (7) applies, and since at optimality problem (12) and (3) have the same objective, we get the result.

Discussion. We can make the following remarks:

1. The per-step cost of ADMM-η is dominated by proxF◦X. For the square loss F(z) = 1 2n ‖z - y‖22 , it requires solving a linear system which is constant for all iterations and hence can be solved efficiently in O(d2 ) time, after a single O(d3 ) operation at the start of the algorithm. Whereas for the hinge loss F(z) = 1n ∑ni=1 max(0, 1 - ziyi), [W. Kienzle, 2006] proposed an "SMO-like" algorithm with O(dn2 ) operations.
2. The assumption make by ADMM-η (regarding proxF◦X) can be compared with those for Alt-η and FISTA-η, both of which requiring proxF◦ ˆX , with ˆX = XD 12 (η). Since the matrix ˆX changes at every iteration, computing proxF◦ ˆX is more difficult than proxF◦X. For instance, for the square loss, one needs to solve a different linear system at each iteration, which amounts to O(d3 ).
3. The step-size selection scheme we described in Section 5.3 leads to the following choices σ = ‖X‖2 /(√3n), τ = √3n/‖X‖2 .

The next section discusses a generic first-order algorithm for (12) when proxF◦X is difficult to compute.

5.3 CP-η: A generic first-order algorithmFor (3)

While proxF◦X may be difficult to obtain, proxF may be easier to compute (this is the case for all common loss functions). We derive a first-order procedure denoted CP-η (Algorithm 3) having a convergence guarantee of O(1/T) for all losses and norms, where proxF and proxΓ are easy to compute.

Algorithm 3 CP-η

Require: X, proxΓ, proxF ∗, k

• 1: Initialize u(0) ∈ R3d , v(0) ∈ Rn+3d (Ref. (20)).
• 2: Initialize r = ‖X‖.
• 3: Define functions G, H∗ as in (21).
• 4: K = [ X 0n×2d -rI3d

]

• 5: (u(k) , v(k) ) =CP(K, proxG, proxH∗, u(0) , v(0) , k)
• 6: w(k) = u(k) (1 : d), η(k) = u(k) (d + 1 : 2d).
• 7: Return w(k) , η(k)

Theorem 2. Let Ψ(w, η) be as defined in (3) with (w∗ , η∗ ) = argminw,η Ψ(w, η), and (u(k) = [w(k) η(k) t(k) ]⊤ , v(k) = [α(k) δ(k) β(k) ν (k) ]⊤ ), K be defined as in Algorithm 3. Then ∀λ > 0, there exists a constant R > 0 such that after T iterations

Proof. (Sketch) Similar to arriving at (6) from (5) through the Fenchel dual of the loss function F, we apply the same trick in (12) to get the following equivalent problem.

For the problem (18), we now use Fenchel duality on the constraint IC and arrive at the following problem.

Now, the mapping of variables in this formulation to Eq. (6) is given by

Here, we have introduced the constant r to balance the scale of α against κ, with the choice r = ‖X‖ in experiments, since ‖K‖2 = ‖X⊤X + r2 I‖2 ⩽ ‖X‖2 + r2 . Algorithm 1 applies and we get Algorithm 3 to solve the problem (12). We see that G is separable function in terms of its variables w, η and t. And so is H∗ separable in terms of α and the triplets (δj, βj, νj) for all j. This makes it easy to compute the proximal operators needed for G and H∗ . Now the result (7) applies, and we get the result.

Discussion. We can make the following remarks:

1. Algorithm 3 is the first ever O(1/T) efficient first-order algorithm for (3), since it accesses X only through matrix-vector multiplications.
2. One may also equate problems (18) and (6) with the following mappings of u, v, G, H ∗ , for which CP achieves O(1/T 2 ) convergence rate when F is smooth (see (8)).

But, proxG is not easy to compute except for simple norms like the ℓ1 or grouped-ℓ1 norm, making the algorithm impractical for general norms.

1. We derive the following step size choices for CP-η: σ = 1/n and τ = n/(2‖X‖2 ).

5.4 General Discussion

We can differentiate the various algorithms for the problem (3) as follows:

• Between Alt-η and FISTA-η, the former does not have any known convergence rates, whereas the latter has a convergence rate of O(1/T 2 ). But the per-step cost of FISTA-η is much higher because of backtracking. As seen in experiments the backtracking cost is quite high which leads to FISTA-η being impractical to use for general losses and norms.
• As discussed before, ADMM-η is similar to FISTA-η and Alt-η, since all of these assume easy computability of proxF◦ ˜X , with ˜X = X for ADMM-η and XD 12 (η) for the other two algorithms. This subtle difference does make an impact, especially for losses like the squared loss, for which the proximal operator for ADMM-η is easier to obtain than those for FISTA-η and Alt-η.
• When proxF◦X is not very efficient to compute, the only choice we have is CP-η, which works with all norms and losses with easy to compute proxF and proxΓ.
• When proxF◦X is easy to compute, the choice of the algorithm depends on the cost of proxΓ. Cheaper the cost of proxΓ computation, better is CP-η in performance compared to ADMM-η. It is because of the experimental observation that ADMM-η takes far fewer iterations than CP-η, and costly proxΓ does make CP-η run longer.
• ADMM-η is observed to be less sensitive to the conditioning of X. This may be because the inner CP routine does not depend on X explicitly, which is taken care of by proxF◦X as a black box. The dependence on X within the CP routine comes only through the step sizes τ and σ. This makes ADMM-η preferred in illconditioned cases.

6 Experiments

To illustrate the efficiency of CP-η and ADMM-η over existing algorithms, we choose the aforementioned tree-sparsity inducing norm ΩH (Example 1), which is popular in wavelet coefficients estimation. As discussed earlier, Ω2S (w) = ΩH(w), but the primal-dual procedures are different because of the difference in the function Γ. Since time complexity to compute proxΓ is identical in both cases (O(d log d)), we choose only ΩH for the purpose of simulations. We include the following solvers: (1) FISTA-η, (2) Alt-η, (3) FISTA 1

1 The non-accelerated variant ISTA was very slow in experiments and hence excluded from the results.

on primal (1) with proxΩ computed using (a) decomposition algorithm (O(d2 )) given in [Obozinski and Bach, 2016], denoted as FISTA-w, (b) Picard iterations [Baldassarre et al., 2012], denoted as FISTA-w(p).

Setup. We perform numerical simulations 2 by generating synthetic data. Following [Bach et al., 2011], we generate X ∈ Rn×d as Xij ∼ N(0, 1). For the ground truth model w∗ , we assumed a tree structure with uniform branching factor of 4 and a depth k. We generated w∗ uniformly at random from [0, 1]d and set s = 0.75d randomly chosen indices to 0 satisfying the tree structure. The labels y were generated as y = Xw + ξ for the squared loss examples, and y = sign(Xw + ξ) for the hinge loss examples, where ξ is a standard Gaussian noise. We fixed n = 1000, d = 15000, λ = 0.01, and the convergence criteria was the relative duality gap (with threshold ϵ = 10-4 ).

6.1 Squared Loss

In this case, proxF◦X can be efficiently computed using matrix-vector multiplications in each iteration, qualifying ADDM-η as a first-order algorithm. We make the following inferences from the simulation plots given in Figure 1.

• All the solutions for the problem (3) do better than that of (1). Both FISTA-w and FISTA-w-p do not converge within a limit of 30 minutes, the latter though being faster than the former suffers from approximate solutions at each step.
• Both CP-η and ADMM-η are better than all the compared algorithms in running time, justifying the claims made in the paper.
• In general ADMM-η is faster than CP-η because of efficient computation of proxF◦X and overall lesser number of iterations.

2 Conducted on a Ubuntu PC with Core i7 processor, 8G RAM.

6.2 Hinge Loss

proxF◦X is computed in this case through solving an SVM [W. Kienzle, 2006], which is based on an SMO algorithm (O(dn2 )). Hence CP-η is the only first-order algorithm in this case. Similar to the squared loss case, both the proposed algorithms perform better than the rest as shown in Figure 2.

6.3 Comparison Of CP-η And ADMM-η

Although Table 2 may suggest that ADMM-η is better than CP-η for both loss functions, this is not always true, as it depends on d and n (see Table 1). When d increases, ADMM-η becomes expensive for square loss due to the O(d2 ) complexity. For hinge loss, both ADMM-η and CP-η have linear dependency on d, and hence theoretically, CP-η has no advantage over ADMM-η in very high dimensional scalings, for a fixed n. However, for ADMM-η, due to quadratic dependency on n (see Table 1), the per-step cost goes high on large sample settings. This dependency on n for ADMMη is justified in Table 2, where we see that the per-step cost for ADMM-η goes much higher for hinge loss compared to square loss.

7 Conclusions and Future Directions

We studied efficient primal-dual algorithms to learn with Subquadratic norms. One may also investigate for potential extensions of these algorithms towards general reweighted least-squares formulations [Bach et al., 2011] for norms. Study of inexact proximal opeartors for subquadratic norms also provides alternate directions.