# Tensorial Model for Photolithography Aerial Image Simulation

Authors: Caroline Fossati; Salah Bourennane; Romuald Sabatier; Antonio Di Giacomo
Journal: Advances in OptoElectronics
 (2009)

---

## Abstract

In this paper, we propose to adapt the multilinear algebra tools to the tensor of Transmission Cross-Coefficients (TCC) values for aerial image simulation in order to keep the data
tensor as a whole entity. This new approach implicitly extends the singular value decomposition (SVD) to tensors, that is, Higher Order SVD or TUCKER3 tensor
decomposition which is used to obtain lower rank-(K1,K2,K3,K4) tensor approximation (LRTA (K1,K2,K3,K4)). This model
requires an Alternating Least Square (ALS) process known as TUCKALS3 algorithm. The needed number of kernels is estimated
using two adapted criteria, well known in signal processing and information theory. For runtime improvement, we use the fixed
point algorithm to calculate only the needed eigenvectors. This new approach leads to a fast and accurate algorithm to compute
aerial images.


---

## Body

## 1. Introduction

Over the last decades, photolithography has been instrumental to the historical trend, better known as Moore’s law, of doubling chip density roughly every eighteen months while maintaining a nearly constant chip price. Going down to the 32 nm node, without changing the wavelength, we move closer and closer to the theoretical optical resolution limit. Therefore resolution Enhancement Techniques (RETs) have been developed in order to period all shape of an integrated circuit properly on the silicon wafer. One of the most well-known RETs is the optical proximity correction (OPC), which consists of making mask pre-compensation of all nonlinear effects, optical diffraction, and interference effects. The reason behind correcting images distortions between the mask and the silicon due to this resolution limit. Model-Based Optical Proximity Correction (MBOPC) can treat layouts by deforming mask pattern to improve resolution on silicon. As for current masks, several billions of segments have to be moved during several iterations necessary to reach convergence; fast and accurate algorithms are mandatory to perform OPC on a mask in a reasonable time for industry.To overcome this limitation, mathematical simplifications have to be done. As imaging with a lithography system (see Figure1) is analogous to microscopy, the theory used in MBOPC is drawn from works originally designed for microscopy theory. Fourier optics was first used to describe the image formed by a microscope [1, 2]. In 1951, Hopkins developed other formulation or method (1), that has been afterward used in other works involving image formation [3–6]. This formulation has the advantage of a four-way transmission function independent of the shape of the object to image. However, this method, in aerial image computation for photolithography, remains highly runtime consuming and needs a numerical approximation in order to reach production constraints.Figure 1
General optical lithography process diagram.In the solution called sum of coherent systems (SOCSs) provided by Cobb in [7], the Transfer Cross-Coefficients (TCCs) (2) are written in matrix format, then Singular Value Decomposition (SVD) that is used on this matrix, resulting in (3), allows to spare a huge amount of time while preserving a good accuracy as only the eigenvectors associated with the most energetic eigenvalues are used in calculus.Clearly, we use Hopkins formulation, and a notation whereX will denote a tensor, X a matrix, x a vector, and x a scalar.The intensity image at the image plane is given by(1)I(u,v)=(12π)2∑i1I1∑i2I2∑i3I3∑i4I4𝒯𝒞𝒞(i1,i2,i3,i4)×E(i1,i2)E*(i3,i4)ei{(i1-i3)u+(i2-i4)v},
where(2)𝒯𝒞𝒞(i1,i2,i3,i4)=2π∬-∞+∞Γ(x,y)F(x+i1,y+i2)×F*(x+i3,y+i4)dxdy,
andI(.) intensity image at the image plane, E(.) object being imaged (mask), Γ(.) mutual intensity function, describing coherence properties of the illumination, F(.) coherent point spread function, describing properties of the projection system.Typically,I1=I2=I3=I4=I for nonastigmatic systems because of symmetry reasons. Those values are spatial frequencies, correspond to the range of values such that 𝒯𝒞𝒞(i1,i2,i3,i4)≠0, and depend on the optical system [8].In his approach, Cobb [7] has unfolded tensor 𝒯𝒞𝒞 to obtain a matrix T representation and then its eigendecomposition (SVD) leads to(3)T=∑k=1I2λkukukH,
where the superscript “H” denotes transpose conjugate. Then the lower rank matrix approximation of T gives the “kernels” used in aerial image calculation.According to (2), the TCCs are depending on four indices, thus in this paper we propose to represent them as a four-way array. In will denote the size of tensor 𝒯𝒞𝒞 in dimension n. Tensor approach permits to treat a four-way array as a whole entity. In the considered problem, this approach allows to preserve interdimensional correlation avoiding possible data loss [9]. Low rank approximation, such as Singular Value Decomposition for two-dimensional array, considers significant and remaining parts of the data, they are based on data of most significant feature selection. The eigendecomposition of the matrix provides eigenvectors [10]. The Kn eigenvectors associated with the Kn largest (dominant) eigenvalues span the so-called “useful subspace.” Useful subspace-based methods are applied to source localization in array processing and image denoising. These methods were adopted to four-way-array (tensor) [9]. The tensor model extends the classical matrix model [9]. This approach implicitly extends the SVD to tensors, that is, Higher Order SVD (HOSVD).A fixed point algorithm comes then replacing HOSVD algorithm to spare some computational load, as already shown in [11].The remainder of the paper is organized as follows In Section2, some tensor definitions and properties are proposed. Section 3 proposes an original lower rank tensor approximation (LRTA) adapted to TCC data approximation. In Section 4, we propose a runtime improvement algorithm using an original method. After presenting some experimental results in Section 5, we conclude the paper in Section 6.


## 2. Tensor Definitions and Properties

In the following we remind some tensor definitions. For more details the reader can be referred to the papers in [12–14].


### 2.1. Tensor Definition

We consider anN-order tensor to represent an N-way array. Each of its elements is accessible through N indices (i1,i2,…,iN). This tensor is written 𝒜∈ℂI1×⋯×IN, and each element ai1,…,iN. Each component is called an “n-mode” referring to the nth index of the tensor. A zero-order tensor is a scalar, a first-order tensor is a vector, and a second-order tensor is a matrix.To ease understanding this paper, some useful definitions are presented in the following.


### 2.2. Unfolding of TensorA

To study multiway data properties in a particularn-mode direction, as shown in Figure 2, we define the n-mode flattening matrix of tensor 𝒜, written An∈ℂIn×Mn, where(4)Mn=I1⋯In-1In+1⋯IN.Figure 2
Flattened matrixA2∈ℝI2×I3I1 in 2-mode of tensor 𝒜∈ℝI1×I2×I3.


### 2.3.n-mode product ×n

Consider tensor𝒜∈ℂI1×⋯×IN and a matrix H∈ℂJ×I with I and J∈ℕ* (natural numbers excluding 0) and for all n∈{1,…,N}, n-mode product between 𝒜 and H is defined as(5)ℬ=𝒜×nH∈ℂI1×⋯In-1×J×In+1⋯×IN,
Where ×n operator defines n-mode product. It generalizes matrix product to tensor and n-mode vector product in a particular n-mode.Each element with index(i1,…,in-1,j,in+1,…,iN) of tensor ℬ is given by(6)bi1…in-1jin+1…iN=∑i=1Inai1…in-1iin+1…iNhji.


### 2.4. Tensorial Scalar Product

Tensorial scalar product is defined as follows. Suppose twoN order tensors: 𝒜 and ℬ∈ℂI1×⋯×IN. Scalar product between 𝒜 and ℬ is given by(7)〈𝒜ℬ〉=∑i1,…,iNai1,…,iNbi1,…,iN.


### 2.5. Useful Subspace Definition

Let us defineE(n) the nth-mode vector space of dimension In, associated with the nth-mode of tensor 𝒜. By definition, E(n) is generated by the column vectors of the nth-mode flattening matrix. The nth-mode flattening matrix An of tensor 𝒜∈ℂI1×⋯×IN is defined as a matrix from ℂIn×Mn, with(8)Mn=I1⋯In-1In+1⋯IN.An columns are the In-dimensional vectors obtained from 𝒜 by varying index in and keeping the other indices fixed. In An, the In-dimensional vectors obtained from 𝒜 are ordered with respect to their index for each mode but the nth mode. These vectors are called in the following the nth-mode vectors of tensor 𝒜. The lower rank-(K1,K2,...,KN) tensor 𝒜 approximation is possible by extending the classical subspace approach to tensors by assuming that, whatever the nth-mode vector space of dimension In, E(n) is the superposition of two orthogonal subspaces:  (i)
the useful subspaceE1(n) of dimension Kn. E1(n) is generated by the Kn eigenvectors of matrix An, which are associated with the Kn largest eigenvalues;(ii)
and the subspaceE2(n) of dimension In-Kn. E2(n) is generated by the In-Kn singular vectors of matrix An, which are associated with the In-Kn smallest eigenvalues,such thatE(n)=E1(n)⊕E2(n). The dimensions K1,K2,…,KN can be estimated by means of the well-known (Akaike Information Criterion) AIC or (Minimum Description Length) MDL criteria [10]. These are entropy-based Information Criteria which have had a fundamental impact in statistical model evaluation problems.


### 2.6. Frobenius Norm of a Tensor

Frobenius norm of a tensor𝒜, written ∥𝒜∥, is given by(9)∥𝒜∥=〈𝒜𝒜〉=∑i1,…,iNai1,…,iN2.


### 2.7. External Product

External product “∘” of several vectors e(1),…,e(N) belonging, respectively, to E(1),…,E(N), of dimension I1,…,IN defines a tensor 𝒜∈ℂI1×…×IN of rank 1:(10)𝒜=e(1)∘e(2)∘⋯∘e(N)
in such a way that whatever the indices (i1,…,iN), the element ai1…iN of tensor 𝒜 is defined by product ai1…iN=∏n=1Nein(n), where ein(n) is the inth component of vector e(n).


### 2.8. Tensor Rank

In case of matrices, rank plays an important role in SVD, in canonical decomposition and in lower rank approximation. Consider the definitions of a tensor rankClassical Rank
Consider tensor of orderN, 𝒜∈ℂI1×⋯×IN. From II-A, tensor 𝒜 is built on N vector space: E(1),…,E(N), of dimension I1,…,IN. Tensor 𝒜 has rank one if, for all n=1,…,N, N vectors u(n)∈E(n) exist, such as 𝒜 is equal to their external product:
(11)𝒜=u(1)∘⋯∘u(N).By extension,𝒜 is a tensor of rank K, if K is the minimal number or tensor of rank 1 which, by summation, give 𝒜. We write the rank of tensor 𝒜 by Rank(𝒜)=K.n-Mode Rank
It is possible to definen-mode rank of a tensor as the generalization of column vectors lines vectors rank of a matrix. n-mode rank of a tensor 𝒜∈ℂI1×⋯×IN, written Rankn(𝒜), is the rank of flattened matrix An of 𝒜 in n-mode, Rankn(𝒜)=Rank(An).Rank-(K1,…,KN)
A tensor𝒜∈ℂI1×⋯×IN has rank-(K1,…,KN) if, for all n=1,…,N, Rankn(𝒜)=Kn.


## 2.1. Tensor Definition

We consider anN-order tensor to represent an N-way array. Each of its elements is accessible through N indices (i1,i2,…,iN). This tensor is written 𝒜∈ℂI1×⋯×IN, and each element ai1,…,iN. Each component is called an “n-mode” referring to the nth index of the tensor. A zero-order tensor is a scalar, a first-order tensor is a vector, and a second-order tensor is a matrix.To ease understanding this paper, some useful definitions are presented in the following.


## 2.2. Unfolding of TensorA

To study multiway data properties in a particularn-mode direction, as shown in Figure 2, we define the n-mode flattening matrix of tensor 𝒜, written An∈ℂIn×Mn, where(4)Mn=I1⋯In-1In+1⋯IN.Figure 2
Flattened matrixA2∈ℝI2×I3I1 in 2-mode of tensor 𝒜∈ℝI1×I2×I3.


## 2.3.n-mode product ×n

Consider tensor𝒜∈ℂI1×⋯×IN and a matrix H∈ℂJ×I with I and J∈ℕ* (natural numbers excluding 0) and for all n∈{1,…,N}, n-mode product between 𝒜 and H is defined as(5)ℬ=𝒜×nH∈ℂI1×⋯In-1×J×In+1⋯×IN,
Where ×n operator defines n-mode product. It generalizes matrix product to tensor and n-mode vector product in a particular n-mode.Each element with index(i1,…,in-1,j,in+1,…,iN) of tensor ℬ is given by(6)bi1…in-1jin+1…iN=∑i=1Inai1…in-1iin+1…iNhji.


## 2.4. Tensorial Scalar Product

Tensorial scalar product is defined as follows. Suppose twoN order tensors: 𝒜 and ℬ∈ℂI1×⋯×IN. Scalar product between 𝒜 and ℬ is given by(7)〈𝒜ℬ〉=∑i1,…,iNai1,…,iNbi1,…,iN.


## 2.5. Useful Subspace Definition

Let us defineE(n) the nth-mode vector space of dimension In, associated with the nth-mode of tensor 𝒜. By definition, E(n) is generated by the column vectors of the nth-mode flattening matrix. The nth-mode flattening matrix An of tensor 𝒜∈ℂI1×⋯×IN is defined as a matrix from ℂIn×Mn, with(8)Mn=I1⋯In-1In+1⋯IN.An columns are the In-dimensional vectors obtained from 𝒜 by varying index in and keeping the other indices fixed. In An, the In-dimensional vectors obtained from 𝒜 are ordered with respect to their index for each mode but the nth mode. These vectors are called in the following the nth-mode vectors of tensor 𝒜. The lower rank-(K1,K2,...,KN) tensor 𝒜 approximation is possible by extending the classical subspace approach to tensors by assuming that, whatever the nth-mode vector space of dimension In, E(n) is the superposition of two orthogonal subspaces:  (i)
the useful subspaceE1(n) of dimension Kn. E1(n) is generated by the Kn eigenvectors of matrix An, which are associated with the Kn largest eigenvalues;(ii)
and the subspaceE2(n) of dimension In-Kn. E2(n) is generated by the In-Kn singular vectors of matrix An, which are associated with the In-Kn smallest eigenvalues,such thatE(n)=E1(n)⊕E2(n). The dimensions K1,K2,…,KN can be estimated by means of the well-known (Akaike Information Criterion) AIC or (Minimum Description Length) MDL criteria [10]. These are entropy-based Information Criteria which have had a fundamental impact in statistical model evaluation problems.


## 2.6. Frobenius Norm of a Tensor

Frobenius norm of a tensor𝒜, written ∥𝒜∥, is given by(9)∥𝒜∥=〈𝒜𝒜〉=∑i1,…,iNai1,…,iN2.


## 2.7. External Product

External product “∘” of several vectors e(1),…,e(N) belonging, respectively, to E(1),…,E(N), of dimension I1,…,IN defines a tensor 𝒜∈ℂI1×…×IN of rank 1:(10)𝒜=e(1)∘e(2)∘⋯∘e(N)
in such a way that whatever the indices (i1,…,iN), the element ai1…iN of tensor 𝒜 is defined by product ai1…iN=∏n=1Nein(n), where ein(n) is the inth component of vector e(n).


## 2.8. Tensor Rank

In case of matrices, rank plays an important role in SVD, in canonical decomposition and in lower rank approximation. Consider the definitions of a tensor rankClassical Rank
Consider tensor of orderN, 𝒜∈ℂI1×⋯×IN. From II-A, tensor 𝒜 is built on N vector space: E(1),…,E(N), of dimension I1,…,IN. Tensor 𝒜 has rank one if, for all n=1,…,N, N vectors u(n)∈E(n) exist, such as 𝒜 is equal to their external product:
(11)𝒜=u(1)∘⋯∘u(N).By extension,𝒜 is a tensor of rank K, if K is the minimal number or tensor of rank 1 which, by summation, give 𝒜. We write the rank of tensor 𝒜 by Rank(𝒜)=K.n-Mode Rank
It is possible to definen-mode rank of a tensor as the generalization of column vectors lines vectors rank of a matrix. n-mode rank of a tensor 𝒜∈ℂI1×⋯×IN, written Rankn(𝒜), is the rank of flattened matrix An of 𝒜 in n-mode, Rankn(𝒜)=Rank(An).Rank-(K1,…,KN)
A tensor𝒜∈ℂI1×⋯×IN has rank-(K1,…,KN) if, for all n=1,…,N, Rankn(𝒜)=Kn.


## 3. Tensor Approximation Based on Lower Rank Tensor Approximation (LRTA) Algorithm

The main objective of this paper is to adopt the multilinear algebra tools to approximate the tensor𝒯𝒞𝒞 for system transmission based on Lower Rank Tensor Approximation (LRTA) algorithm; we will consider in the rest of the paper N=4 which corresponds to the fourth-order tensor of the TCC values, 𝒯𝒞𝒞∈ℂI1×I2×I3×I4. The corresponding n-mode rank values are denoted K1,K2,K3,K4. We assume that the value of Kn is known or is estimated (see Section 4.2 to know how to estimate them) for all n=1 to 4.As an extension to the vector and matrix cases, in the tensor formulation, the projectors on thenth-mode vector spaces are determined by computing the rank-(K1,K2,K3,K4) approximation of 𝒯𝒞𝒞 in the least-squares sense. From a mathematical point of view, the rank-(K1,K2,K3,K4) approximation of 𝒯𝒞𝒞 is represented by tensor 𝒯𝒞𝒞̂ which minimizes the quadratic tensor Frobenius norm ∥𝒯𝒞𝒞-𝒯𝒞𝒞̂∥2 subject to the constraint that 𝒯𝒞𝒞̂∈ℂI1×I2×I3×I4 is a rank-(K1,K2,K3,K4) tensor. The TUCKER3 tensor decomposition also known as Higher-Order SVD (HOSVD) [15] and lower rank-(K1,K2,K3,K4) tensor approximation (LRTA-(K1,K2,K3,K4)) has recently been used as multimode (Principal Component Analysis) PCA, in seismic for wave separation based on a useful subspace method, in image processing for face recognition and expression analysis and noise filtering of color images. Tensor 𝒯𝒞𝒞̂ can be expressed with the TUCKER3 model [12] as(12)𝒯𝒞𝒞̂=𝒟×1U(1)×2U(2)×3U(3)×4U(4).
The least square solution implies that(13)𝒟=𝒯𝒞𝒞×1U(1)T×2U(2)T×3U(3)T×4U(4)T.Thus the best lower approximation of tensor𝒯𝒞𝒞 denoted 𝒯𝒞𝒞̂ is to orthogonally project, for every n-mode, the vectors of tensor 𝒯𝒞𝒞 on the nth-mode useful subspace E1(n), for all n=1 to 4 [9, 13](14)𝒯𝒞𝒞̂=𝒯𝒞𝒞×1P(1)×2P(2)×3P(3)×4P(4).In (14), P(n)=U(n)U(n)T is called the projector of the nth-mode space. The estimation of the projectors is a difficult nonlinear least square problem that is generally solved owing to an ALS algorithm referred to as TUCKALS3 algorithm [12, 16].TUCKALS3 algorithm [12, 16] is an optimal algorithm to estimate the different n-mode useful subspaces. The description of TUCKALS3 algorithm, used in rank-(K1,K2,K3,K4) approximation, is provided in Algorithm 1.Algorithm 1:Lower rank-(K1,K2,K3,K4) tensor approximation (LRTA-(K1,K2,K3,K4))
(1)Input: data tensor 𝒯𝒞𝒞, and dimensions K1,K2,K3,K4 of all n-mode signal subspaces.(2)Initialization  k=0: For n=1 to 4, calculate the projectors P0(n) given by HOSVD-(K1,K2,K3,K4):(a)n-mode unfold 𝒯𝒞𝒞 into matrix TCCn;(b) Compute the SVD ofTCCn;(c) Compute matrixU0(n) formed by the Kn eigenvectors associated with the Kn largest singular values ofTCCn. U0(n) is the initial matrix of the n-mode signal subspace orthogonal basis vectors;(d) Form the initial orthogonal projectorP0(n)=U0(n)U0(n)T on the n-mode signal subspace;(e) Compute the HOSVD-(K1,K2,K3,K4) of tensor 𝒯𝒞𝒞 given byℬ0=𝒯𝒞𝒞×1P0(1)×2⋯×4P0(4);(3)ALS loop: Repeat until convergence, that is, for example, while ∥ℬk+1-ℬk∥2>ε, ε>0 being a priorfixed threshold,(a) Forn=1 to 4:(i) Formℬ(n),k: ℬ(n),k=𝒯𝒞𝒞×1Pk+1(1)×2Pk+1(2)×3Pk+1(3)×4Pk+1(4);(ii)n-mode unfold tensor ℬ(n),k into matrix Bn(n),k;(iii) Compute matrixC(n),k=Bn(n),kTCCnT;(iv) Compute matrixUk+1(n) composed of the Kn eigenvectors associated with the Kn largest eigenvaluesofC(n),k. Uk(n) is the matrix of the n-mode signal subspace orthogonal basis vectors at the kthiteration;(v) ComputePk+1(n)=Uk+1(n)Uk+1(n)T;(b) Computeℬk+1=𝒯𝒞𝒞×1Pk+1(1)×2⋯×4Pk+1(4);(c) Incrementk.(4)Output: the estimated signal tensor is obtained through 𝒯𝒞𝒞̂=𝒯𝒞𝒞×1Pkstop(1)×2⋯×4Pkstop(4). 𝒯𝒞𝒞̂ isthe rank-(K1,K2,K3,K4) approximation of 𝒯𝒞𝒞, where kstop is the index of the last iteration after theconvergence of TUCKALS3 algorithm.In this algorithm, the second-order statistics comes from the SVD of matrixTCCn at step 2b, which is equivalent, up to 1/Mn multiplicative factor, to the estimation of tensor 𝒯𝒞𝒞n-mode vectors [9]. The definition of Mn is given in (8). In the same way, at step 3(a)iii, matrix C(n),k is, up to 1/Mn multiplicative factor, the estimation of the covariance matrix between tensor 𝒯𝒞𝒞 and tensor ℬ(n),kn-mode vectors. According to step 3(a)i, ℬ(n),k represents data tensor 𝒯𝒞𝒞 filtered in every mth-mode but the nth-mode, by projection-filters Pl(m), with m≠n, l=k if m>n and l=k+1 if m<n. TUCKALS3 algorithm has recently been used to process a multimode PCA in order to perform white noise removal in color images [9]. A good approximation of the rank-(K1,K2,K3,K4) approximation can simply be achieved by computing the  HOSVD-(K1,K2,K3,K4) of tensor 𝒯𝒞𝒞 [9, 15]. Indeed, the HOSVD-(K1,K2,K3,K4) of 𝒯𝒞𝒞 consists of the initialization step of TUCKALS3 algorithm, and hence can be considered as a suboptimal solution for the rank-(K1,K2,K3,K4) approximation of tensor 𝒯𝒞𝒞 [15]. This HOSVD-based technique has recently been used in [9] for image denoising.


## 4. Runtime Improvement: Fixed Point Algorithm

As computing LRTA involves to compute four times Singular Value Decomposition for𝒯𝒞𝒞 decomposition, we propose a fixed point algorithm based on Gram-Schmidt orthogonalisation to lower runtime.


### 4.1. Fixed Point Algorithm

We propose to adapt a Fixed Point algorithm [17] to replace SVD in LRTA [11], Algorithm 1, it allows to spare a high amount of computational load while preserving a high level of accuracy and has never been used in optical domain yet. Although SVD needs to be computed only once, this method is very useful during model tuning as several variants are compiled.One way to compute theK orthonormal basis vectors is to use Gram-Schmidt method (Algorithm 2).Algorithm 2:Fixed-point
(1)  ChooseK, the number of principal axes or eigenvectors required to estimate. Consider matrix T and setp←1.(2)  Initialize eigenvectorup of size d×1, for example, randomly;(3)  Updateup as up←Tup;(4)  Do the Gram-Schmidt orthogonalization processup←up-∑j=1j=p-1(upTuj)uj;(5)  Normalizeup by dividing it by its norm: up←up/∥up∥.(6)  Ifup has not converged, go back to step 3.(7)  Increment counterp←p+1 and go to step 2 until p equals K.The eigenvector with dominant eigenvalue will be measured first. Similarly, all the remainingK-1 basis vectors (orthonormal to the previously measured basis vectors) will be measured one by one in a reducing order of dominance. The previously measured (p-1)th basis vectors will be used to find the pth basis vector. The algorithm for pth basis vector will converge when the new value up+ and old value up are such that up+Tup=1. It is usually economical to use a finite tolerance error to satisfy the convergence criterion ∥up+Tup-1∥<η, where η is a prior fixed threshold.LetU=[u1,u2,…,uK] be the matrix whose columns are the K orthonormal basis vectors. Then UUT is the projector onto the subspace spanned by the K eigenvectors associated with the K largest eigenvalues. This subspace is also called “signal subspace.” It can be used in LRTA-(K1,K2,K3,K4) to retrieve the basis vectors U0(n) in step 2c of Algorithm 1. Thus, the initialization step is faster since it does not need the In basis vectors but only the Kn first ones and it does not need the step 2b, that is SVD of the data tensor n-mode unfolding matrix TCCn.According to what have been explained before, fixed point algorithm will replace SVD in LRTA algorithm, in order to improve runtime. However, a prerequisite to this algorithm is the knowledge of eachn-mode rank. A method to determine them is developed in next subsection.


### 4.2.n-Mode Rank Estimation

Whereas a  common way to obtainK1,K2,K3,K4 values is to use empirical data [18], we wish to present here a scientific tool optimizing the calculation. Each projector P(n), n=1,2,3,4, is estimated by the truncation of unfolding matrices TCCn using SVD, that is, by keeping the Kn eigenvectors associated with the Kn largest singular values of TCCn,n=1,2,3,4. In order to estimate the value of Kn for each n-mode, we extend the well-known detection criteria [10]. These criteria, initially, are developed for estimating the number of largest eigenvalues of any matrix. Thus, the optimal signal subspace dimension is obtained merely by minimizing one of (Akaike’s Information Criterion) AIC [19] or (Minimum Description Length) MDL [20] criteria. This approach does note require any subjective threshold settings. The number of eigenvectors can be determined from the eigenvalues of the n-mode unfolding matrix TCCn. Thus, this method presents the advantage of being self-consistent as it computes the optimal number of kernels for any illumination configuration.Consequently, for eachn-mode unfolding of 𝒯𝒞𝒞, the form of detection criterion AIC can be expressed as(15)AIC(k)=-2(In-k)Mnln𝒢(λk+1,…,λIn)𝒜(λk+1,…,λIn)+2k(2In-k),
and the MDL criterion is given by(16)MDL(k)=-(In-k)Mnln𝒢(λk+1,…,λIn)𝒜(λk+1,…,λIn)+12k(2In-k)lnMn,
where 𝒢 and 𝒜 are geometricand arithmetic mean, respectively. (λi)1≤i≤In are In eigenvalues of the data matrix of the n-mode unfolding 𝒯𝒞𝒞: λ1≥λ2≥⋯≥λIn, and Mn is the number of columns of the n-mode unfolding 𝒯𝒞𝒞.Then-mode rank Kn is the value of k (k=1,…,In-1) which minimizes AIC or MDL criterion.


## 4.1. Fixed Point Algorithm

We propose to adapt a Fixed Point algorithm [17] to replace SVD in LRTA [11], Algorithm 1, it allows to spare a high amount of computational load while preserving a high level of accuracy and has never been used in optical domain yet. Although SVD needs to be computed only once, this method is very useful during model tuning as several variants are compiled.One way to compute theK orthonormal basis vectors is to use Gram-Schmidt method (Algorithm 2).Algorithm 2:Fixed-point
(1)  ChooseK, the number of principal axes or eigenvectors required to estimate. Consider matrix T and setp←1.(2)  Initialize eigenvectorup of size d×1, for example, randomly;(3)  Updateup as up←Tup;(4)  Do the Gram-Schmidt orthogonalization processup←up-∑j=1j=p-1(upTuj)uj;(5)  Normalizeup by dividing it by its norm: up←up/∥up∥.(6)  Ifup has not converged, go back to step 3.(7)  Increment counterp←p+1 and go to step 2 until p equals K.The eigenvector with dominant eigenvalue will be measured first. Similarly, all the remainingK-1 basis vectors (orthonormal to the previously measured basis vectors) will be measured one by one in a reducing order of dominance. The previously measured (p-1)th basis vectors will be used to find the pth basis vector. The algorithm for pth basis vector will converge when the new value up+ and old value up are such that up+Tup=1. It is usually economical to use a finite tolerance error to satisfy the convergence criterion ∥up+Tup-1∥<η, where η is a prior fixed threshold.LetU=[u1,u2,…,uK] be the matrix whose columns are the K orthonormal basis vectors. Then UUT is the projector onto the subspace spanned by the K eigenvectors associated with the K largest eigenvalues. This subspace is also called “signal subspace.” It can be used in LRTA-(K1,K2,K3,K4) to retrieve the basis vectors U0(n) in step 2c of Algorithm 1. Thus, the initialization step is faster since it does not need the In basis vectors but only the Kn first ones and it does not need the step 2b, that is SVD of the data tensor n-mode unfolding matrix TCCn.According to what have been explained before, fixed point algorithm will replace SVD in LRTA algorithm, in order to improve runtime. However, a prerequisite to this algorithm is the knowledge of eachn-mode rank. A method to determine them is developed in next subsection.


## 4.2.n-Mode Rank Estimation

Whereas a  common way to obtainK1,K2,K3,K4 values is to use empirical data [18], we wish to present here a scientific tool optimizing the calculation. Each projector P(n), n=1,2,3,4, is estimated by the truncation of unfolding matrices TCCn using SVD, that is, by keeping the Kn eigenvectors associated with the Kn largest singular values of TCCn,n=1,2,3,4. In order to estimate the value of Kn for each n-mode, we extend the well-known detection criteria [10]. These criteria, initially, are developed for estimating the number of largest eigenvalues of any matrix. Thus, the optimal signal subspace dimension is obtained merely by minimizing one of (Akaike’s Information Criterion) AIC [19] or (Minimum Description Length) MDL [20] criteria. This approach does note require any subjective threshold settings. The number of eigenvectors can be determined from the eigenvalues of the n-mode unfolding matrix TCCn. Thus, this method presents the advantage of being self-consistent as it computes the optimal number of kernels for any illumination configuration.Consequently, for eachn-mode unfolding of 𝒯𝒞𝒞, the form of detection criterion AIC can be expressed as(15)AIC(k)=-2(In-k)Mnln𝒢(λk+1,…,λIn)𝒜(λk+1,…,λIn)+2k(2In-k),
and the MDL criterion is given by(16)MDL(k)=-(In-k)Mnln𝒢(λk+1,…,λIn)𝒜(λk+1,…,λIn)+12k(2In-k)lnMn,
where 𝒢 and 𝒜 are geometricand arithmetic mean, respectively. (λi)1≤i≤In are In eigenvalues of the data matrix of the n-mode unfolding 𝒯𝒞𝒞: λ1≥λ2≥⋯≥λIn, and Mn is the number of columns of the n-mode unfolding 𝒯𝒞𝒞.Then-mode rank Kn is the value of k (k=1,…,In-1) which minimizes AIC or MDL criterion.


## 5. Experimental Results

The proposed methods can be applied to any four-way array of data set such as image, multicomponent seismic signals, hyperspectral images.We wish to compare LRTA algorithm benefits over SOCS algorithm. We use square size matrices with increasing rows and columns number as data set to compare algorithms runtime and accuracy.


### 5.1. Runtime Improvement

We first focus on Fixed Point algorithm runtime gain toward SVD algorithm. Figure3 shows runtime curves of increasing size matrices decomposition. Matlab SVD algorithm is used. We have to compute about 10 eigenvectors for each matrix with fixed point algorithm. This value is provided by AIC criterion and corresponds to an average value of eigenvectors number used in classical OPC models. However, as fixed point algorithm runtime is linear with eigenvectors number to compute, Table 1 shows that this number can be increased up to around 300 eigenvectors for 3000×3000 matrices and up to around 100 eigenvectors for 120×120 matrices.Table 1
Some numerical values corresponding with Figure3.
Matrix size1207201320192025203000SVD runtime (sec)0.1096.26238.911122.810279.762489.808Fixed point runtime (sec)0.0100.4061.7624.5369.22014.522SVD/FP rate10.62915.44122.08327.07330.34333.728Figure 3
Computational times (in sec.) as a function of rows and columns number.The computations have been run with Matlab on a 2.4 GHz dual core Pentium with 4Go RAM under Windows XP.


### 5.2. AIC and MDL Criteria

Figure4(a) shows AIC and MDL functions of a fourth-order tensor filled with random data. The minimum value in AIC and MDL curves gives rank estimation value of considered mode. Figure 4(b) shows the normalized energy of the eigenvalues of the same array. Those curves are helpful to understand that the rank given by AIC and MDL criteria is really connected to the eigenvectors associated with the largest eigenvalues.Energy of matrix eigenvalues and comparison between AIC and MDL criteria.
(a)
AIC and MDL estimation criteria(b)
Normalized energy of matrix eigenvalues


### 5.3. Fixed Point Algorithm Reconstruction Error

In order to assess our algorithm precision toward SVD decomposition, Figure5 shows reconstruction error of both fixed point and SVD algorithms, ∥T-T′∥/∥T∥, where T denotes the original matrix and T′ the reconstructed matrix. For constancy purposes with fixed point algorithm, the same number of eigenvectors is used for both algorithms. It is difficult to distinguish the two curves on Figure 5; however fixed point algorithm reconstruction error remains very slightly above SVD reconstruction error.Figure 5
Reconstruction error comparison between SVD and Fixed point algorithm.


### 5.4. Comparison between LRTA and SOCS Algorithms

Figure6 shows a comparison of SOCS and LRTA algorithms reconstruction error and both algorithms runtime.Comparison between LRTA and SOCS.
(a)
Comparison of LRTA and SOCS reconstruction error(b)
Comparison between SOCS and LRTA runtimeTo obtain these curves, random tensors of increasing size have been created. They have been unfolded following Cobb’s method [7], decomposed using SVD algorithm, recomposed using only the dominant eigenvectors.In the same time, they have been decomposed using LRTA algorithm with fixed point and also recomposed using only the dominant eigenvectors.In both cases, the error function is the same than used in Section5.3 to compute fixed point algorithm versus SVD algorithm error.The curves show clearly the benefit of our method based on tensor data representation and on multilinear algebra tools.The computations have been run withMatlab on a 2.4 GHz dual core Pentium with 4Go RAM under Windows XP.


### 5.5. Comparison of Aerial Images

Figure7 shows aerial images computed using Full Hopkins equation and with SOCS algorithm and LRTA algorithm, for 4 (Figure 7(a)), 10 (Figure 7(b)) and 20 kernels (Figure 7(c)). The structure used is a T-shaped geometry of total size 720nm×720nm. The model used corresponds to a 193 nm wavelength illumination, 0.75 NA and 0.6 σ, optical diameter is 1280 nm, and grid size is 1 nm.Aerial image comparison between SOCS, LRTA, and full Hopkins.
(a)(b)(c)A lower number of kernels are sufficient for the LRTA algorithm to give a result close to the full hopkins. This is confirmed by the obtained results shown in Figure8 using the same calculus only with 4 kernels for L-shaped geometry. Note that this sufficient number is estimated using AIC criterion.Figure 8
Aerial image comparison between SOCS, LRTA, and full Hopkins.Considering Full Hopkins gives the aerial image as close as possible to the image really obtained on the wafer. One can see that the aerial image obtained using LRTA approximation better looks like real image than one obtained by SOCS algorithm.


## 5.1. Runtime Improvement

We first focus on Fixed Point algorithm runtime gain toward SVD algorithm. Figure3 shows runtime curves of increasing size matrices decomposition. Matlab SVD algorithm is used. We have to compute about 10 eigenvectors for each matrix with fixed point algorithm. This value is provided by AIC criterion and corresponds to an average value of eigenvectors number used in classical OPC models. However, as fixed point algorithm runtime is linear with eigenvectors number to compute, Table 1 shows that this number can be increased up to around 300 eigenvectors for 3000×3000 matrices and up to around 100 eigenvectors for 120×120 matrices.Table 1
Some numerical values corresponding with Figure3.
Matrix size1207201320192025203000SVD runtime (sec)0.1096.26238.911122.810279.762489.808Fixed point runtime (sec)0.0100.4061.7624.5369.22014.522SVD/FP rate10.62915.44122.08327.07330.34333.728Figure 3
Computational times (in sec.) as a function of rows and columns number.The computations have been run with Matlab on a 2.4 GHz dual core Pentium with 4Go RAM under Windows XP.


## 5.2. AIC and MDL Criteria

Figure4(a) shows AIC and MDL functions of a fourth-order tensor filled with random data. The minimum value in AIC and MDL curves gives rank estimation value of considered mode. Figure 4(b) shows the normalized energy of the eigenvalues of the same array. Those curves are helpful to understand that the rank given by AIC and MDL criteria is really connected to the eigenvectors associated with the largest eigenvalues.Energy of matrix eigenvalues and comparison between AIC and MDL criteria.
(a)
AIC and MDL estimation criteria(b)
Normalized energy of matrix eigenvalues


## 5.3. Fixed Point Algorithm Reconstruction Error

In order to assess our algorithm precision toward SVD decomposition, Figure5 shows reconstruction error of both fixed point and SVD algorithms, ∥T-T′∥/∥T∥, where T denotes the original matrix and T′ the reconstructed matrix. For constancy purposes with fixed point algorithm, the same number of eigenvectors is used for both algorithms. It is difficult to distinguish the two curves on Figure 5; however fixed point algorithm reconstruction error remains very slightly above SVD reconstruction error.Figure 5
Reconstruction error comparison between SVD and Fixed point algorithm.


## 5.4. Comparison between LRTA and SOCS Algorithms

Figure6 shows a comparison of SOCS and LRTA algorithms reconstruction error and both algorithms runtime.Comparison between LRTA and SOCS.
(a)
Comparison of LRTA and SOCS reconstruction error(b)
Comparison between SOCS and LRTA runtimeTo obtain these curves, random tensors of increasing size have been created. They have been unfolded following Cobb’s method [7], decomposed using SVD algorithm, recomposed using only the dominant eigenvectors.In the same time, they have been decomposed using LRTA algorithm with fixed point and also recomposed using only the dominant eigenvectors.In both cases, the error function is the same than used in Section5.3 to compute fixed point algorithm versus SVD algorithm error.The curves show clearly the benefit of our method based on tensor data representation and on multilinear algebra tools.The computations have been run withMatlab on a 2.4 GHz dual core Pentium with 4Go RAM under Windows XP.


## 5.5. Comparison of Aerial Images

Figure7 shows aerial images computed using Full Hopkins equation and with SOCS algorithm and LRTA algorithm, for 4 (Figure 7(a)), 10 (Figure 7(b)) and 20 kernels (Figure 7(c)). The structure used is a T-shaped geometry of total size 720nm×720nm. The model used corresponds to a 193 nm wavelength illumination, 0.75 NA and 0.6 σ, optical diameter is 1280 nm, and grid size is 1 nm.Aerial image comparison between SOCS, LRTA, and full Hopkins.
(a)(b)(c)A lower number of kernels are sufficient for the LRTA algorithm to give a result close to the full hopkins. This is confirmed by the obtained results shown in Figure8 using the same calculus only with 4 kernels for L-shaped geometry. Note that this sufficient number is estimated using AIC criterion.Figure 8
Aerial image comparison between SOCS, LRTA, and full Hopkins.Considering Full Hopkins gives the aerial image as close as possible to the image really obtained on the wafer. One can see that the aerial image obtained using LRTA approximation better looks like real image than one obtained by SOCS algorithm.


## 6. Conclusion

We have proposed here a new approach to Transfer Cross-Coefficients data set approximation in diffraction theory of optical images based on multilinear algebra tools. The goal of this work is two way axed, it allows to adapt complex physical equations to tensor computation and to improve runtime by using fixed point algorithm instead of HOSVD while preserving accuracy by using LRTA and AIC algorithm. We have proven our method to be faster and more accurate than existing SOCS method, as we obtain a lower reconstruction error than SOCS algorithm, whatever the size of input data is. Those methods have proven their efficiency in other domains, such as four-way images restoration or denoising and for wave separation of multicomponent seismic data, but have never been implemented in imaging theory yet. In further work, we will bring more detailed experimental data to express in real cases the extent of the improvement of our method.



---

*Source: 457549-2009-12-06.xml*