# Graph-Based Symbolic Technique and Its Application in the Frequency Response Bound Analysis of Analog Integrated Circuits

Authors: E. Tlelo-Cuautle; S. Rodriguez-Chavez; A. A. Palma-Rodriguez
Journal: The Scientific World Journal
 (2014)

---

## Abstract

A new graph-based symbolic technique (GBST) for deriving exact analytical expressions like the transfer functionH(s) of an analog integrated circuit (IC), is introduced herein. The derived H(s) of a given analog IC is used to compute the frequency response bounds (maximum and minimum) associated to the magnitude and phase of H(s), subject to some ranges of process variational parameters, and by performing nonlinear constrained optimization. Our simulations demonstrate the usefulness of the new GBST for deriving the exact symbolic expression for H(s), and the last section highlights the good agreement between the frequency response bounds computed by our variational analysis approach versus traditional Monte Carlo simulations. As a conclusion, performing variational analysis using our proposed GBST for computing the frequency response bounds of analog ICs, shows a gain in computing time of 100x for a differential circuit topology and 50x for a 3-stage amplifier, compared to traditional Monte Carlo simulations.


---

## Body

## 1. Introduction

With the downscaling of the integrated circuit (IC) technology, nanometer circuit designs become more and more sensitive to process variations [1, 2], which are produced by fluctuations at the moment of manufacturing, and have been continuously increasing in relative magnitude as IC technology continues to scale to 45 nm and below. On the one hand, IC designers usually perform repeated Monte Carlo (MC) simulations to predict variations, leading to an expensive computational cost, while the main goal is computing the bounds of a given performance by varying the value of some parameters under certain percentage. On the other hand, to mitigate this drawback on performing repeated MC simulations, performance bound methods emerged as attractive techniques for variational analysis of analog ICs under parameter variations [3–7]. However, although those bound methods are quite efficient computing the lower and upper bounds of the magnitude and phase of a transfer function H(s), systematic methods have not been proposed yet to obtain the variational symbolic expression for H(s). Henceforth, this paper introduces a performance bound technique that derives the exact analytical H(s) of a linear (linearized) analog IC, as CMOS amplifiers, by applying a new graph-based symbolic technique (GBST), which is described in the next section. Afterwards, the frequency response bounds, for the magnitude and phase for H(s) of three amplifiers, are obtained by applying nonlinear constrained optimization.


## 2. Graph-Based Symbolic Technique

Symbolic analysis of analog circuits has been demonstrated to be useful for deriving analytical expressions such asH(s) [8–14]. An important contribution on the development of symbolic analysis techniques was the introduction of determinant decision diagrams (DDDs) [8, 9], whose advantage is computing symbolic expressions sharing many subexpressions [10, 11]. However, it was not realized for deriving exact analytical expressions for CMOS amplifiers [15]. That way, we introduce herein a new graph-based symbolic technique (GBST) for deriving exact symbolic expressions of analog ICs. Our proposed GBST exploits DDD concepts and has the advantage of being compact and unique as well.


### 2.1. Simple Case: Symbolic Determinant without Node Reuse

Let us consider the determinant given by (1) [9], of size n×n. By applying GBST, the graph representation is built in a depth-first search (DFS) fashion, while one expects having paths of n+1 levels. Every element in the graph corresponds to a nonzero entry in M. In this manner, one obtains the graph shown in Figure 1, where applying the rule of signs from Cramer’s rule (see (2)) does the assignation of signs to each node. A path is eliminated if a zero entry in the nodal admittance matrix is found. Consider
(1)|M|=|ab00cde00fgh00ij|=adgj-adhi-aefj-bcgj+bchi,(2)sign⁡=(-1)row+col.Figure 1
Graph representation of (1).A tree in which arithmetic operations are encoded in the depth of the tree nodes represents the graph. That is, different depth implies multiplication and equal depth implies addition. This leads us to derive the expression:(3)|M|=a[d(gj-hi)+e(-fj)]+b[c(-gj+hi)].


### 2.2.  Advanced Case: Symbolic Determinant with Node Reuse

In (3), the determinant is expressed as a sequence of products [8, 9], because the graph shown in Figure 1 does not reuse node information, while having many repetitive terms that correspond to the repeated minors gj and hi. For this case, the smallest the matrix minor, the highest the repetition rate. Therefore, our main idea is reusing the information of those repeated nodes. For example: gj, hi, and fj are product terms of 2×2 minors in (1). In this manner, there are five two-nodded subgraphs with vertex sets V1={g,j}, V2={h,i}, V3={f,j}, V4={g,j}, and V5={h,i}. As one sees, edges E1={g,j} and E4={g,j} are equivalents, so that subgraphs G1={V1,E1} and G4={V4,E4} carry the same information. In this case, node reuse is possible, and one just needs to identify node ancestors; for example, the subgraphs G1=G4 and G3={V3,E3} have the node ancestors: g for G1=G4 and f for G3. Recall that each node is linked to a nonzero matrix entry; that is, Arow,col≠0. Node f is in turn a representation of A3,2 and g of A3,3.Extending node reuse in the whole graph, the first obvious consequence is that there are no repeated nodes; in other words, for a matrixA with nz nonzero entries, there are nz nodes.


### 2.3.  The Advanced Case in Detail

To derive the symbolic expression with node reuse, three different data structures are required. The first and most obvious is the node structure that contains the following:(i)
node name: a unique name for each node, assigned as an index number;(ii)
terms: an array containing the index and sign of the element;(iii)
column: the column of the nonzero entry where the node belongs to;(iv)
descendants: an array of node pointers linking to the descendants of the current node in the graph structure.The second data structure is a graph type with the fields:(i)
graph name: a unique name for the current graph. It is possible to have many different graphs; for example, to compute the transfer functionH(s), two graphs are required: one for deriving the numerator and the other for the denominator [9];(ii)
matrix size: the size of the square matrix;(iii)
root node: it is a trivial node with term value equal to 1 (see Figure1) and row and column equal to zero. When multiple graphs are constructed (during factorization), the root node can be any of the nz nodes;(iv)
visited columns: when traversing a graph to represent the determinant, the column of a visited node is appended to this array.The third structure stores matrix elements and the independent vector.From these structures,nz+1 nodes are created. The row and column fields are used to evaluate the sign by (2) as well as to determine which nodes are to be skipped. The graph is generated starting with a trivial node named 0 with term value of 1; then the multiplication of terms is codified as the depth in the graph, and nodes are linked accordingly. The algorithm to build the graph structure for the representation of |A| is sketched by Algorithm 1. Algorithm 2 computes the expression for the determinant. That way, the graph associated to the determinant in (1) is shown in Figure 2.Algorithm 1:Build graph (A(i,j), Ancestor).
m
← number of Columns of A
n
← number of Rows of A
D
← set of descendants of Ancestor
Ensure: m>0∧n=m
function buildGraph(A, Ancestor)
if  m>1  then
prune=1
for  i=1 to m  do
for all Columns j that Ai,j≠0  do
if BUILDGRAPH(Cij,Ai,j) ≠0  then
D=D∪Ai,j
prune=0
else
if  Ai,j≠0  then
prune=0
return  pruneAlgorithm 2:Symbolic determinant from a graph.
(1) Read the graph in a DFS fashion:
(a) Keep track of which columns have been visited
(b) Skip nodes from columns already visited
(2) The symbolic expression is the sum of products according to the visited descendantsFigure 2
Determinant graph with node reuse.


### 2.4. Symbol Factorization

Factorization takes place by executing Algorithm3. The result is the expression as a polynomial that is represented by an array of sum of products with one entry for each power of symbol w.Algorithm 3:Graph to polynomial form.
(1) Expand the nodes in the graph containing the symbolw.
(2) Read the graph in a DFS fashion and preserve only those routes from
root to bottom with at least one occurrence of symbolw.
(3) The number of occurrences of the symbol is the power ofw for a given route.
(4) Express each route as a sum of products.
(5) The coefficients are the summation of all appended routes for each power ofw.Let us consider matrixM in (4). To derive the polynomial form as powers of w, the associated graph is shown in Figure 3, and the polynomial is given by (5). Consider
(4)M=[a+wb00cde00fgh+w00ij],(5)w0=adgj-adhi-aefj-bcgj+bchi,w1=adi+dgj-dhi-efj-bci,w2=di.Figure 3
Expanded graph.


### 2.5. Applying the Graph-Based Symbolic Technique

Our proposed GBST starts off with an SPICE-like netlist as input. The allowed circuit elements are R, C, L, V, I, E, G, and M, being resistor, capacitor, inductor, independent voltage source, independent current source, voltage-controlled voltage source, voltage-controlled current source, and MOSFET. When formulating the equations, the symbol name is taken exactly as specified in the netlist, for example, R_name, C_name, M_name, and so forth. Their corresponding small-signal model based on controlled sources substitutes all active elements. In turn, controlled sources are modeled with combinations of norator and nullator in order to make use of the extensive studies on analysis of Nullor based circuits given in [9, 11, 14, 15]. The guidelines for obtaining the nodal admittance matrix by applying nodal analysis are summarized in [11]. Afterwards, applying our proposed graph-based symbolic technique derives the solution. In this manner, the computing information on the solutions for the common source amplifier, differential pair amplifier, three stages uncompensated OTA [15], recycled folded cascode OTA [16], and 741 opamp is listed in Table 1. Using the values computed by HSPICE, in order to verify correctness, does the numerical evaluation of the derived symbolic expression.Table 1
Symbolic formulation and numerical evaluation ofD(s),N(s),andH(s) by applying the new GBST.
Circuit features
Computer time (seconds)
Circuit
Elements
Nodes
Equations
D(s)
N(s)
H(s)
Differential pair
35
26
1.1235
0.122
0.1464
1.4895
RFC OTA [16]
106
56
1.6603
0.201
0.1869
2.2633
LV Amp [15]
33
18
2.35
0.058
0.0464
2.4544
Common source
8
6
0.8581
0.041
0.0205
0.9811
741
112
77
0.5123
1.37
0.822
2.7043To demonstrate the suitability of the new graph-based symbolic technique (GBST), we performed a comparison with HSPICE and the DDD symbolic tool [3], as shown in Figure 4. As one sees, GBST is in good agreement with the numerical response computed by SPICE, while the DDD technique has an error around 10%. That way, GBST is applied herein to derive the exact analytical expression H(s) of CMOS amplifiers. The derived H(s) is used to perform variational analysis, in order to compute the frequency response bounds (maximum and minimum) of the magnitude and phase from H(s), subject to some ranges of process variational parameters by performing nonlinear constrained optimization.Figure 4
H
(
s
) comparison among DDD [3] (dashed), HSPICE (dots), and our graph-based symbolic technique (solid) for the differential pair topology.


## 2.1. Simple Case: Symbolic Determinant without Node Reuse

Let us consider the determinant given by (1) [9], of size n×n. By applying GBST, the graph representation is built in a depth-first search (DFS) fashion, while one expects having paths of n+1 levels. Every element in the graph corresponds to a nonzero entry in M. In this manner, one obtains the graph shown in Figure 1, where applying the rule of signs from Cramer’s rule (see (2)) does the assignation of signs to each node. A path is eliminated if a zero entry in the nodal admittance matrix is found. Consider
(1)|M|=|ab00cde00fgh00ij|=adgj-adhi-aefj-bcgj+bchi,(2)sign⁡=(-1)row+col.Figure 1
Graph representation of (1).A tree in which arithmetic operations are encoded in the depth of the tree nodes represents the graph. That is, different depth implies multiplication and equal depth implies addition. This leads us to derive the expression:(3)|M|=a[d(gj-hi)+e(-fj)]+b[c(-gj+hi)].


## 2.2.  Advanced Case: Symbolic Determinant with Node Reuse

In (3), the determinant is expressed as a sequence of products [8, 9], because the graph shown in Figure 1 does not reuse node information, while having many repetitive terms that correspond to the repeated minors gj and hi. For this case, the smallest the matrix minor, the highest the repetition rate. Therefore, our main idea is reusing the information of those repeated nodes. For example: gj, hi, and fj are product terms of 2×2 minors in (1). In this manner, there are five two-nodded subgraphs with vertex sets V1={g,j}, V2={h,i}, V3={f,j}, V4={g,j}, and V5={h,i}. As one sees, edges E1={g,j} and E4={g,j} are equivalents, so that subgraphs G1={V1,E1} and G4={V4,E4} carry the same information. In this case, node reuse is possible, and one just needs to identify node ancestors; for example, the subgraphs G1=G4 and G3={V3,E3} have the node ancestors: g for G1=G4 and f for G3. Recall that each node is linked to a nonzero matrix entry; that is, Arow,col≠0. Node f is in turn a representation of A3,2 and g of A3,3.Extending node reuse in the whole graph, the first obvious consequence is that there are no repeated nodes; in other words, for a matrixA with nz nonzero entries, there are nz nodes.


## 2.3.  The Advanced Case in Detail

To derive the symbolic expression with node reuse, three different data structures are required. The first and most obvious is the node structure that contains the following:(i)
node name: a unique name for each node, assigned as an index number;(ii)
terms: an array containing the index and sign of the element;(iii)
column: the column of the nonzero entry where the node belongs to;(iv)
descendants: an array of node pointers linking to the descendants of the current node in the graph structure.The second data structure is a graph type with the fields:(i)
graph name: a unique name for the current graph. It is possible to have many different graphs; for example, to compute the transfer functionH(s), two graphs are required: one for deriving the numerator and the other for the denominator [9];(ii)
matrix size: the size of the square matrix;(iii)
root node: it is a trivial node with term value equal to 1 (see Figure1) and row and column equal to zero. When multiple graphs are constructed (during factorization), the root node can be any of the nz nodes;(iv)
visited columns: when traversing a graph to represent the determinant, the column of a visited node is appended to this array.The third structure stores matrix elements and the independent vector.From these structures,nz+1 nodes are created. The row and column fields are used to evaluate the sign by (2) as well as to determine which nodes are to be skipped. The graph is generated starting with a trivial node named 0 with term value of 1; then the multiplication of terms is codified as the depth in the graph, and nodes are linked accordingly. The algorithm to build the graph structure for the representation of |A| is sketched by Algorithm 1. Algorithm 2 computes the expression for the determinant. That way, the graph associated to the determinant in (1) is shown in Figure 2.Algorithm 1:Build graph (A(i,j), Ancestor).
m
← number of Columns of A
n
← number of Rows of A
D
← set of descendants of Ancestor
Ensure: m>0∧n=m
function buildGraph(A, Ancestor)
if  m>1  then
prune=1
for  i=1 to m  do
for all Columns j that Ai,j≠0  do
if BUILDGRAPH(Cij,Ai,j) ≠0  then
D=D∪Ai,j
prune=0
else
if  Ai,j≠0  then
prune=0
return  pruneAlgorithm 2:Symbolic determinant from a graph.
(1) Read the graph in a DFS fashion:
(a) Keep track of which columns have been visited
(b) Skip nodes from columns already visited
(2) The symbolic expression is the sum of products according to the visited descendantsFigure 2
Determinant graph with node reuse.


## 2.4. Symbol Factorization

Factorization takes place by executing Algorithm3. The result is the expression as a polynomial that is represented by an array of sum of products with one entry for each power of symbol w.Algorithm 3:Graph to polynomial form.
(1) Expand the nodes in the graph containing the symbolw.
(2) Read the graph in a DFS fashion and preserve only those routes from
root to bottom with at least one occurrence of symbolw.
(3) The number of occurrences of the symbol is the power ofw for a given route.
(4) Express each route as a sum of products.
(5) The coefficients are the summation of all appended routes for each power ofw.Let us consider matrixM in (4). To derive the polynomial form as powers of w, the associated graph is shown in Figure 3, and the polynomial is given by (5). Consider
(4)M=[a+wb00cde00fgh+w00ij],(5)w0=adgj-adhi-aefj-bcgj+bchi,w1=adi+dgj-dhi-efj-bci,w2=di.Figure 3
Expanded graph.


## 2.5. Applying the Graph-Based Symbolic Technique

Our proposed GBST starts off with an SPICE-like netlist as input. The allowed circuit elements are R, C, L, V, I, E, G, and M, being resistor, capacitor, inductor, independent voltage source, independent current source, voltage-controlled voltage source, voltage-controlled current source, and MOSFET. When formulating the equations, the symbol name is taken exactly as specified in the netlist, for example, R_name, C_name, M_name, and so forth. Their corresponding small-signal model based on controlled sources substitutes all active elements. In turn, controlled sources are modeled with combinations of norator and nullator in order to make use of the extensive studies on analysis of Nullor based circuits given in [9, 11, 14, 15]. The guidelines for obtaining the nodal admittance matrix by applying nodal analysis are summarized in [11]. Afterwards, applying our proposed graph-based symbolic technique derives the solution. In this manner, the computing information on the solutions for the common source amplifier, differential pair amplifier, three stages uncompensated OTA [15], recycled folded cascode OTA [16], and 741 opamp is listed in Table 1. Using the values computed by HSPICE, in order to verify correctness, does the numerical evaluation of the derived symbolic expression.Table 1
Symbolic formulation and numerical evaluation ofD(s),N(s),andH(s) by applying the new GBST.
Circuit features
Computer time (seconds)
Circuit
Elements
Nodes
Equations
D(s)
N(s)
H(s)
Differential pair
35
26
1.1235
0.122
0.1464
1.4895
RFC OTA [16]
106
56
1.6603
0.201
0.1869
2.2633
LV Amp [15]
33
18
2.35
0.058
0.0464
2.4544
Common source
8
6
0.8581
0.041
0.0205
0.9811
741
112
77
0.5123
1.37
0.822
2.7043To demonstrate the suitability of the new graph-based symbolic technique (GBST), we performed a comparison with HSPICE and the DDD symbolic tool [3], as shown in Figure 4. As one sees, GBST is in good agreement with the numerical response computed by SPICE, while the DDD technique has an error around 10%. That way, GBST is applied herein to derive the exact analytical expression H(s) of CMOS amplifiers. The derived H(s) is used to perform variational analysis, in order to compute the frequency response bounds (maximum and minimum) of the magnitude and phase from H(s), subject to some ranges of process variational parameters by performing nonlinear constrained optimization.Figure 4
H
(
s
) comparison among DDD [3] (dashed), HSPICE (dots), and our graph-based symbolic technique (solid) for the differential pair topology.


## 3. Constrained Optimization

To highlight the appropriateness of applying the new GBST in variational analysis for finding the lower and upper performance bounds for the magnitude and phase of an amplifier, we formulate a nonlinear constrained optimization problem. That way, we start from a transfer function, whose coefficients (e.g.,ai(p1,…,pm) and bj(p1,…,pn)) are obtained by our proposed graph-based symbolic technique in ans-expanded form; that is,
(6)H(ω)=∑i=0mai(p1,…,pm)si∑j=0nbj(p1,…,pn)sj.Notice thatH(s,p1,…,pm) is a nonlinear function of p1,…,pm. Furthermore, each parameter pi is a random variable with a variational range.Fors=jω, the nominal transfer function becomes
(7)H(s)=H(jω)=H(jω)ejθ(ω).From it, one can obtain a variational transfer function with bounded magnitude and phase regions, which is described by [3]
(8)Hl(ω)≤H(ω)≤Hu(ω),θl(ω)≤θ(ω)≤θu(ω),whereHl(ω) and Hu(ω) are the lower and upper bounds of the magnitude, respectively, and θl(ω) and θu(ω) are the lower and upper bounds of the phase. The evaluation of (7) gives a complex valued result, where the magnitude H(ω)=|H(jω)| and the phase angle θ(ω)=∠H(jω) are real values. The goal is to find the bounds of the magnitude and phase for H(jω), such that one can obtain (8). Henceforth, in the presence of process variations, the signal is perturbed from its nominal behavior, and it is usually bounded between its minimum and maximum limits, as sketched in Figure 5.Figure 5
Transfer function description including parameter variations.For instance, using the lower bound of the magnitude response|H(jω)| at frequency ω, then the objective function is minimized and described by
(9)minimize|H(jω,x)|subjecttoxl≤x≤xu.In (9) x=[p1,…,pm] represents the circuit parameters variable vector, subject to the optimization constraints [xl,xu]. In IC design foundries and cell library vendors supply the constraints. Algorithm 4 summarizes this process.Algorithm 4:Calculation of frequency response bounds using our graph-based symbolic technique from Section 2 and constrained optimization.
(1) Parse circuit netlist.
(2) Set bounds on parameters affected by process variations.
(3) Derive the symbolic expression, for example, the transfer function.
(4)for All frequencies ωi  do
(5)    Constrained optimization to find magnitude bounds onωi.


### 3.1. A Simple Example in MATLAB

Let us consider Figure6. In the frequency domain, capacitors and inductors are analyzed as complex impedances, and then the transfer function is given by
(10)H(jω)=VoutVin=1/jωCR+jωL+(1/jωC).Figure 6
RLC circuit.Assuming thatC and L have variations of 20% from their nominal values C=1μF and L=1μH, then C∈[0.8,1.2]μF and L∈[0.8,1.2]μH. For this example, the iterative method called active set was used. As a result, three snapshots at different frequency points are shown in Figure 7.Figure 7
Frequency response of the RLC circuit [3]. Solid curve is the magnitude response with nominal parameters. The dashed curves are the lower and upper bounds due to variations. The three surfaces at the top have L and C as x-axis and y-axis, respectively, and z-axis shows the magnitude variations at three sampling frequencies.


### 3.2. Line Search Algorithm

A line search algorithm searches for a decrease inf in a descent direction using the Armijo rule [17, 18] for stepsize control. The steplength is chosen to minimize f along the ray {x+λΔx∣λ≥0}, where Δx is called thedescent direction and the point x+λΔx is called x+, or even to just reduce f “enough” [19]. That way, given the current point xc and descent direction Δx, it looks for λ, such that
(11)f(xc+λΔx)≤f(xc).However, if the decreasing achieved by this inequality for someλ is too small, it is not possible to guarantee convergence to a local minimum. So in order to avoid this issue, λ must satisfyArmijo rule (sufficient decrease) given by
(12)f(xc+λΔx)≤f(xc)+λγ∇f(xc)TΔx,
where γ∈(0,1) [20]. This is shown in Figure 8.Figure 8
Line search representation.In Figure8, the interval where (12) is accomplished is [0,λmax⁡]. We can rewrite this by (13) where the Armijo condition [18] to accept a trial point x+ is given by (14), and M is an integer greater than zero. If x+ is rejected, the steplength is redefined by (15), where λ* minimizes a quadratic model. This strategy of repeatedly testing sufficient decrease and reducing the stepsize if the test fails is called backtracking. The projected gradient (PGRAD) and spectral projected gradient (SPG) methods used herein for variational analysis of CMOS amplifiers are based on backtracking line search. Consider
(13)f(x+)≤f(xc)+γ∇f(xc)T(x+-xc),(14)f(x+)≤max⁡0≤j≤min⁡⁡{c,M-1}f(xc-j)+γ∇f(xc)T(x+-xc),(15)λ=max⁡{σ1λ,min⁡{σ2λ,λ*}}.


### 3.3. Projected Gradient Method

The PGRAD method is described in [21]. It is the natural extension of thesteepest descent algorithm [20], used in unconstrained optimization to bound constrained problems.Given a current iteratexc, the new iterate is
(16)x+=P(xc-λ∇f(xc)),
where the gradient ∇f(x) is defined by
(17)∇f(x)=(∂f(x)∂x1,…,∂f(x)∂xn)hhhhhhforx=(x1,…,xn).In (17), λ is a steplength parameter given by the Armijo rule, where one searches on a ray from xc in a direction where f is locally decreasing [18]. Besides, in order to implement a line search scheme, one must specify what sufficient decrease means. Therefore, for bound constrained problems thesufficient decrease condition for line search is defined by
(18)f(x(λ))-f(x)≤-αλ∥x-x(λ)∥2,
where
(19)x(λ)=P(x-λ∇f(x))forλ≥0.In this caseλ=βm, where β∈(0,1) and m≥0, which is the smallest nonnegative integer such that there issufficient decrease in f. For the end condition, it is necessary to define the active and inactive sets. The set of constraints is calledfeasible set (Ω), and a point in this set is afeasible point. Because thefeasible set is compact, there is always a solution for this minimization problem. The ith constraint isactive at x∈Ω if either xi=Li or xi=Ui. If the ith constraint is notactive, it is calledinactive. Therefore, one can write A(x) and I(x) for the active and inactive sets, respectively. Anactive/inactive set is the set of indexes i such that the ith constraint is active/inactive. In this manner,
(20)r0=∥x0-x0(1)∥.The stop criterion is given by (21), where τr and τa are relative and absolute tolerances. Consider
(21)∥x-x(1)∥≤τa+τrr0.This optimization method is shown in Algorithm5. It starts with x0∈Ω and uses a sufficient decrease parameter γ∈(0,1) and safeguarding parameters 0≤σ1≤σ2≤1. Initially, α0=1/∥f(x0)∥2. Given xc and αc≥0 the algorithm shows how to obtain xc+1 and αk+1 and when to stop the optimization engine.Algorithm 5:Projected gradient method.
(1)if  ∥P(xc-∇f(xc))-xc∥=0, that is, xc is stationary then
(2)   Stop
(3) Computedc=P(xc-αc∇f(xc))-xc
(4) Setλ←1
(5)  Setx+=x++λxc
(6)if  f(x+)≤f(xc)+γ〈x+-xc,∇f(xc)〉  then
(7)λc=λ, xc+1=x+, sk=xc+1-xc, yk=∇f(xc+1)
(8)   Go to Step 12
(9)else
(10)  Defineλnew∈[σ1λ,σ2λ] and λ←λnew
(11)  Go to Step 5
(12) Compute the trial step lengthαc+1=αc∥f(xc)∥2/∥f(xc+1)∥2


### 3.4. Spectral Projected Gradient Method

The SPG method improves choosing the steplength, which greatly speeds up the convergence of gradient methods [22]. Unlike the projected gradient method, the spectral projected gradient one [21] is more related to the quasi-Newton family [17]. The main idea behind the spectral choice of steplength is that the steepest descent method is very slow but it can be accelerated taking, instead of the stepsize that comes from the minimization of the function along the gradient of the current iteration, the one that comes from the one-dimensional minimization at the previous step.The point in the first iteration of this method should be a feasible point; that is, the algorithm starts withx0∈Ω and uses an integer M≥1, a small parameter αmin⁡≥0, a large parameter αmax⁡≥αmin⁡, a sufficient decrease parameter γ∈(0,1), and safeguarding parameters 0≤σ1≤σ2≤1. Consider ∥P(xc-∇f(xc))-xc∥=0.The method is shown in Algorithm6, where λnew uses one-dimensional quadratic interpolation and it is safeguarded taking λ←λ/2 when the minimum of the one-dimensional quadratic lies outside [0.1,0.9λ]. The line search condition in step 6 guarantees that the sequence {xc} remains in Ω0={x∈Ω:f(x)≤f(x0)}.Algorithm 6:Spectral projected gradient method.
(1)if  ∥P(xc-∇f(xc))-xc∥=0, that is, xc is stationary then
(2)   Stop
(3) Computedc=P(xc-αcg(xc))-xc.
(4) Setλ←1
(5)   Setx+=xc+λΔf(xc)
(6)if  f(x+)≤max⁡0≤j≤min⁡{c,M-1}f(xc-j)+γλ〈dc,∇f(xc)〉  then
(7)   Defineλc=λ, xc+1=x+, sc=xc+1-xc, yc=∇f(xc+1)-∇f(xc)
(8)   Got to Step 11
(9)else
(10)  Setλnew∈[σ1λ,σ2λ], λ←λnew  and go to Step 5
(11)  Computebc=〈sc,yc〉.
(12)if  bk≤0  then
(13)αc+1=αmax⁡.
(14)else
(15)αc=〈sc,sc〉 and αc+1=min⁡{αmax⁡,max⁡{αmin⁡,αc/bc}}


### 3.5. Remarks

The projected gradient and spectral projected gradient methods start atx0∈Ω and use as search direction the internal projected gradient direction. In case of rejection of the first trial point, the next ones are computed along the same line. Also, for both methods, the calculation of λnew uses a one-dimensional quadratic interpolation. Both algorithms involve a projection on the convex set Ω, a function evaluation f(x), and a gradient evaluation per iteration ∇f(x).


## 3.1. A Simple Example in MATLAB

Let us consider Figure6. In the frequency domain, capacitors and inductors are analyzed as complex impedances, and then the transfer function is given by
(10)H(jω)=VoutVin=1/jωCR+jωL+(1/jωC).Figure 6
RLC circuit.Assuming thatC and L have variations of 20% from their nominal values C=1μF and L=1μH, then C∈[0.8,1.2]μF and L∈[0.8,1.2]μH. For this example, the iterative method called active set was used. As a result, three snapshots at different frequency points are shown in Figure 7.Figure 7
Frequency response of the RLC circuit [3]. Solid curve is the magnitude response with nominal parameters. The dashed curves are the lower and upper bounds due to variations. The three surfaces at the top have L and C as x-axis and y-axis, respectively, and z-axis shows the magnitude variations at three sampling frequencies.


## 3.2. Line Search Algorithm

A line search algorithm searches for a decrease inf in a descent direction using the Armijo rule [17, 18] for stepsize control. The steplength is chosen to minimize f along the ray {x+λΔx∣λ≥0}, where Δx is called thedescent direction and the point x+λΔx is called x+, or even to just reduce f “enough” [19]. That way, given the current point xc and descent direction Δx, it looks for λ, such that
(11)f(xc+λΔx)≤f(xc).However, if the decreasing achieved by this inequality for someλ is too small, it is not possible to guarantee convergence to a local minimum. So in order to avoid this issue, λ must satisfyArmijo rule (sufficient decrease) given by
(12)f(xc+λΔx)≤f(xc)+λγ∇f(xc)TΔx,
where γ∈(0,1) [20]. This is shown in Figure 8.Figure 8
Line search representation.In Figure8, the interval where (12) is accomplished is [0,λmax⁡]. We can rewrite this by (13) where the Armijo condition [18] to accept a trial point x+ is given by (14), and M is an integer greater than zero. If x+ is rejected, the steplength is redefined by (15), where λ* minimizes a quadratic model. This strategy of repeatedly testing sufficient decrease and reducing the stepsize if the test fails is called backtracking. The projected gradient (PGRAD) and spectral projected gradient (SPG) methods used herein for variational analysis of CMOS amplifiers are based on backtracking line search. Consider
(13)f(x+)≤f(xc)+γ∇f(xc)T(x+-xc),(14)f(x+)≤max⁡0≤j≤min⁡⁡{c,M-1}f(xc-j)+γ∇f(xc)T(x+-xc),(15)λ=max⁡{σ1λ,min⁡{σ2λ,λ*}}.


## 3.3. Projected Gradient Method

The PGRAD method is described in [21]. It is the natural extension of thesteepest descent algorithm [20], used in unconstrained optimization to bound constrained problems.Given a current iteratexc, the new iterate is
(16)x+=P(xc-λ∇f(xc)),
where the gradient ∇f(x) is defined by
(17)∇f(x)=(∂f(x)∂x1,…,∂f(x)∂xn)hhhhhhforx=(x1,…,xn).In (17), λ is a steplength parameter given by the Armijo rule, where one searches on a ray from xc in a direction where f is locally decreasing [18]. Besides, in order to implement a line search scheme, one must specify what sufficient decrease means. Therefore, for bound constrained problems thesufficient decrease condition for line search is defined by
(18)f(x(λ))-f(x)≤-αλ∥x-x(λ)∥2,
where
(19)x(λ)=P(x-λ∇f(x))forλ≥0.In this caseλ=βm, where β∈(0,1) and m≥0, which is the smallest nonnegative integer such that there issufficient decrease in f. For the end condition, it is necessary to define the active and inactive sets. The set of constraints is calledfeasible set (Ω), and a point in this set is afeasible point. Because thefeasible set is compact, there is always a solution for this minimization problem. The ith constraint isactive at x∈Ω if either xi=Li or xi=Ui. If the ith constraint is notactive, it is calledinactive. Therefore, one can write A(x) and I(x) for the active and inactive sets, respectively. Anactive/inactive set is the set of indexes i such that the ith constraint is active/inactive. In this manner,
(20)r0=∥x0-x0(1)∥.The stop criterion is given by (21), where τr and τa are relative and absolute tolerances. Consider
(21)∥x-x(1)∥≤τa+τrr0.This optimization method is shown in Algorithm5. It starts with x0∈Ω and uses a sufficient decrease parameter γ∈(0,1) and safeguarding parameters 0≤σ1≤σ2≤1. Initially, α0=1/∥f(x0)∥2. Given xc and αc≥0 the algorithm shows how to obtain xc+1 and αk+1 and when to stop the optimization engine.Algorithm 5:Projected gradient method.
(1)if  ∥P(xc-∇f(xc))-xc∥=0, that is, xc is stationary then
(2)   Stop
(3) Computedc=P(xc-αc∇f(xc))-xc
(4) Setλ←1
(5)  Setx+=x++λxc
(6)if  f(x+)≤f(xc)+γ〈x+-xc,∇f(xc)〉  then
(7)λc=λ, xc+1=x+, sk=xc+1-xc, yk=∇f(xc+1)
(8)   Go to Step 12
(9)else
(10)  Defineλnew∈[σ1λ,σ2λ] and λ←λnew
(11)  Go to Step 5
(12) Compute the trial step lengthαc+1=αc∥f(xc)∥2/∥f(xc+1)∥2


## 3.4. Spectral Projected Gradient Method

The SPG method improves choosing the steplength, which greatly speeds up the convergence of gradient methods [22]. Unlike the projected gradient method, the spectral projected gradient one [21] is more related to the quasi-Newton family [17]. The main idea behind the spectral choice of steplength is that the steepest descent method is very slow but it can be accelerated taking, instead of the stepsize that comes from the minimization of the function along the gradient of the current iteration, the one that comes from the one-dimensional minimization at the previous step.The point in the first iteration of this method should be a feasible point; that is, the algorithm starts withx0∈Ω and uses an integer M≥1, a small parameter αmin⁡≥0, a large parameter αmax⁡≥αmin⁡, a sufficient decrease parameter γ∈(0,1), and safeguarding parameters 0≤σ1≤σ2≤1. Consider ∥P(xc-∇f(xc))-xc∥=0.The method is shown in Algorithm6, where λnew uses one-dimensional quadratic interpolation and it is safeguarded taking λ←λ/2 when the minimum of the one-dimensional quadratic lies outside [0.1,0.9λ]. The line search condition in step 6 guarantees that the sequence {xc} remains in Ω0={x∈Ω:f(x)≤f(x0)}.Algorithm 6:Spectral projected gradient method.
(1)if  ∥P(xc-∇f(xc))-xc∥=0, that is, xc is stationary then
(2)   Stop
(3) Computedc=P(xc-αcg(xc))-xc.
(4) Setλ←1
(5)   Setx+=xc+λΔf(xc)
(6)if  f(x+)≤max⁡0≤j≤min⁡{c,M-1}f(xc-j)+γλ〈dc,∇f(xc)〉  then
(7)   Defineλc=λ, xc+1=x+, sc=xc+1-xc, yc=∇f(xc+1)-∇f(xc)
(8)   Got to Step 11
(9)else
(10)  Setλnew∈[σ1λ,σ2λ], λ←λnew  and go to Step 5
(11)  Computebc=〈sc,yc〉.
(12)if  bk≤0  then
(13)αc+1=αmax⁡.
(14)else
(15)αc=〈sc,sc〉 and αc+1=min⁡{αmax⁡,max⁡{αmin⁡,αc/bc}}


## 3.5. Remarks

The projected gradient and spectral projected gradient methods start atx0∈Ω and use as search direction the internal projected gradient direction. In case of rejection of the first trial point, the next ones are computed along the same line. Also, for both methods, the calculation of λnew uses a one-dimensional quadratic interpolation. Both algorithms involve a projection on the convex set Ω, a function evaluation f(x), and a gradient evaluation per iteration ∇f(x).


## 4. Simulation Results

The analog ICs for testing the variational methods PGRAD and SPG that useH(s) derived by the new GBST are shown in Figure 9. That way, after deriving the transfer function by applying the new GBST, and after setting the desired frequencies, the PGRAD and SPG methods are applied to find the magnitude and phase bounds for H(s).Circuits used for the variational analysis: (a) differential pair, (b) 2-stage OTA, and (c) 3-stage OTA.
(a)
(b)
(c)Both the PGRAD and SPG methods were programmed in C and compiled in an Ubuntu Linux environment with the GNU C compiler gcc-4.6.1, with 4 GB RAM in an Intel Core i3. The performance bounds results for the differential pair are shown in Figure10 and for the 3-stage OTA in Figure 11. The lines in blue are those corresponding to Monte Carlo simulations, which are well bounded by the bounds computed by the variational methods PGRAD and SPG. Table 2 summarizes the results where it can be appreciated that both variational methods show better times than by performing repeated Monte Carlo simulations when using HSPICE.Table 2
Time comparison between the variational methods and HSPICE.
Circuit
PGRAD
SPG
HSPICE
Differential pair
15.312 ms
16.363 ms
1.716 s
2-stage OTA
44.559 ms
91.936 ms
81.65 ms
3-stage OTA
38.315 ms
39.233 ms
2.163 sBounds for the differential pair applying (a) PGRAD and (b) SPG.
(a)
(b)Bounds for the 3-stage OTA applying (a) PGRAD and (b) SPG.
(a)
(b)


## 5. Conclusions

A new graph-based symbolic technique (GBST) for deriving exact analytical expressions for analog ICs, such asH(s), and whose symbolic expressions are used to perform variational analysis to obtain performance frequency bounds has been introduced. The variational analysis needs as input the H(s) computed by applying our proposed GBST, and it was based on nonlinear constrained optimization using two line search methods, namely, projected gradient (PGRAD) and spectral projected gradient (SPG).It was demonstrated that the bound analysis based on GBST is suitable for analog ICs. In addition, the time computation of PGRAD and SPG for obtaining the frequency response bounds of analog ICs was improved with respect to HSPICE, when performing repeated Monte Carlo simulations. As shown in the last section, computing the frequency response bounds by using GBST showed a gain of 100x for the differential pair and 50x for the 3-stage amplifier, compared to repeated Monte Carlo simulations, thus justifying the usefulness of using GBST.



---

*Source: 202371-2014-07-17.xml*