Tensorial Model for Photolithography Aerial Image Simulation

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-(𝐾₁,𝐾₂,𝐾₃,𝐾₄) tensor approximation (LRTA (𝐾₁,𝐾₂,𝐾₃,𝐾₄)). 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.

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 Figure 1) 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 where 𝑋 will denote a tensor, 𝑋 a matrix, 𝑥 a vector, and 𝑥 a scalar.

The intensity image at the image plane is given by1𝐼(𝑢,𝑣)=()2𝜋^2𝐼1_𝑖1^𝐼2_𝑖2^𝐼3_𝑖3^𝐼4_𝑖4𝒯𝒞𝒞(𝑖₁,𝑖₂,𝑖₃,𝑖₄)×𝐸(𝑖₁,𝑖₂)𝐸^∗(𝑖₃,𝑖₄)𝑒^{𝑖{(𝑖1−𝑖3)𝑢+(𝑖2−𝑖4)𝑣}}, (1) where𝒯𝒞𝒞(𝑖₁,𝑖₂,𝑖₃,𝑖₄)=2𝜋^+∞_−∞Γ(𝑥,𝑦)𝐹(𝑥+𝑖₁,𝑦+𝑖₂)×𝐹^∗(𝑥+𝑖₃,𝑦+𝑖₄)𝑑𝑥𝑑𝑦, (2) and

𝐼(.) intensity image at the image plane, 𝐸(.) object being imaged (mask), Γ(.) mutual intensity function, describing coherence properties of the illumination, 𝐹(.) coherent point spread function, describing properties of the projection system.

Typically, 𝐼₁=𝐼₂=𝐼₃=𝐼₄=𝐼 for nonastigmatic systems because of symmetry reasons. Those values are spatial frequencies, correspond to the range of values such that 𝒯𝒞𝒞(𝑖₁,𝑖₂,𝑖₃,𝑖₄)≠0, and depend on the optical system [8].

In his approach, Cobb [7] has unfolded tensor 𝒯𝒞𝒞 to obtain a matrix 𝐓 representation and then its eigendecomposition (SVD) leads to𝐓=^𝐼2_𝑘=1𝜆_𝑘𝐮_𝑘𝐮^𝐻_𝑘, (3) where the superscript “H” denotes transpose conjugate. Then the lower rank matrix approximation of 𝐓 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. 𝐼_𝑛 will denote the size of tensor 𝒯𝒞𝒞 in dimension 𝑛. 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 𝐾_𝑛 eigenvectors associated with the 𝐾_𝑛 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 Section 2, 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 an 𝑁-order tensor to represent an N-way array. Each of its elements is accessible through 𝑁 indices (𝑖₁,𝑖₂,…,𝑖_𝑁). This tensor is written 𝒜∈ℂ^{𝐼1×⋯×𝐼𝑁}, and each element 𝑎_{𝑖1,…,𝑖𝑁}. Each component is called an “𝑛-mode” referring to the 𝑛th 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 Tensor 𝐴

To study multiway data properties in a particular 𝑛-mode direction, as shown in Figure 2, we define the 𝑛-mode flattening matrix of tensor 𝒜, written 𝐀_𝑛∈ℂ^{𝐼𝑛×𝑀𝑛}, where𝑀_𝑛=𝐼₁⋯𝐼_𝑛−1𝐼_𝑛+1⋯𝐼_𝑁. (4)

Figure 2 

Flattened matrix 𝐀2∈ℝ𝐼2×𝐼3𝐼1 in 2-mode of tensor 𝒜∈ℝ𝐼1×𝐼2×𝐼3.

2.3. 𝑛-mode product ×_𝑛

Consider tensor 𝒜∈ℂ^{𝐼1×⋯×𝐼𝑁} and a matrix 𝐇∈ℂ^𝐽×𝐼 with 𝐼 and 𝐽∈ℕ^∗ (natural numbers excluding 0) and for all 𝑛∈{1,…,𝑁}, 𝑛-mode product between 𝒜 and 𝐇 is defined asℬ=𝒜×_𝑛𝐇∈ℂ^{𝐼1×⋯𝐼𝑛−1×𝐽×𝐼𝑛+1⋯×𝐼𝑁}, (5) Where ×_𝑛 operator defines 𝑛-mode product. It generalizes matrix product to tensor and 𝑛-mode vector product in a particular 𝑛-mode.

Each element with index (𝑖₁,…,𝑖_𝑛−1,𝑗,𝑖_𝑛+1,…,𝑖_𝑁) of tensor ℬ is given by𝑏_{𝑖1…𝑖𝑛−1𝑗𝑖𝑛+1…𝑖𝑁}=^𝐼𝑛_𝑖=1𝑎_{𝑖1…𝑖𝑛−1𝑖𝑖𝑛+1…𝑖𝑁}ℎ_𝑗𝑖. (6)

2.4. Tensorial Scalar Product

Tensorial scalar product is defined as follows. Suppose two 𝑁 order tensors: 𝒜 and ℬ∈ℂ^{𝐼1×⋯×𝐼𝑁}. Scalar product between 𝒜 and ℬ is given by⟨𝒜ℬ⟩=_{𝑖1,…,𝑖𝑁}𝑎_{𝑖1,…,𝑖𝑁}𝑏_{𝑖1,…,𝑖𝑁}. (7)

2.5. Useful Subspace Definition

Let us define 𝐸^(𝑛) the 𝑛th-mode vector space of dimension 𝐼_𝑛, associated with the 𝑛th-mode of tensor 𝒜. By definition, 𝐸^(𝑛) is generated by the column vectors of the 𝑛th-mode flattening matrix. The 𝑛th-mode flattening matrix 𝐀_𝑛 of tensor 𝒜∈ℂ^{𝐼1×⋯×𝐼𝑁} is defined as a matrix from ℂ^{𝐼𝑛×𝑀𝑛}, with𝑀_𝑛=𝐼₁⋯𝐼_𝑛−1𝐼_𝑛+1⋯𝐼_𝑁. (8)

𝐀_𝑛 columns are the 𝐼_𝑛-dimensional vectors obtained from 𝒜 by varying index 𝑖_𝑛 and keeping the other indices fixed. In 𝐀_𝑛, the 𝐼_𝑛-dimensional vectors obtained from 𝒜 are ordered with respect to their index for each mode but the 𝑛th mode. These vectors are called in the following the 𝑛th-mode vectors of tensor 𝒜. The lower rank-(𝐾₁,𝐾₂,...,𝐾_𝑁) tensor 𝒜 approximation is possible by extending the classical subspace approach to tensors by assuming that, whatever the 𝑛th-mode vector space of dimension 𝐼_𝑛, 𝐸^(𝑛) is the superposition of two orthogonal subspaces: (i) the useful subspace 𝐸₁^(𝑛) of dimension 𝐾_𝑛. 𝐸₁^(𝑛) is generated by the 𝐾_𝑛 eigenvectors of matrix 𝐀_𝑛, which are associated with the 𝐾_𝑛 largest eigenvalues; (ii) and the subspace 𝐸₂^(𝑛) of dimension 𝐼_𝑛−𝐾_𝑛. 𝐸₂^(𝑛) is generated by the 𝐼_𝑛−𝐾_𝑛 singular vectors of matrix 𝐀_𝑛, which are associated with the 𝐼_𝑛−𝐾_𝑛 smallest eigenvalues,

such that 𝐸^(𝑛)=𝐸₁^(𝑛)⊕𝐸₂^(𝑛). The dimensions 𝐾₁,𝐾₂,…,𝐾_𝑁 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‖𝒜‖=⟨𝒜𝒜⟩=_{𝑖1,…,𝑖𝑁}𝑎²_{𝑖1,…,𝑖𝑁}. (9)

2.7. External Product

External product “∘” of several vectors 𝐞⁽¹⁾,…,𝐞^(𝑁) belonging, respectively, to 𝐸⁽¹⁾,…,𝐸^(𝑁), of dimension 𝐼₁,…,𝐼_𝑁 defines a tensor 𝒜∈ℂ^{𝐼1×…×𝐼𝑁} of rank 1:𝒜=𝐞⁽¹⁾∘𝐞⁽²⁾∘⋯∘𝐞^(𝑁) (10) in such a way that whatever the indices (𝑖₁,…,𝑖_𝑁), the element 𝑎_{𝑖1…𝑖𝑁} of tensor 𝒜 is defined by product 𝑎_{𝑖1…𝑖𝑁}=∏^𝑁_𝑛=1𝑒_𝑖^(𝑛)_𝑛, where 𝑒_𝑖^(𝑛)_𝑛 is the 𝑖^th_𝑛 component of vector 𝐞^(𝑛).

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 rank

Classical Rank
Consider tensor of order 𝑁, 𝒜∈ℂ^{𝐼1×⋯×𝐼𝑁}. From II-A, tensor 𝒜 is built on 𝑁 vector space: 𝐸⁽¹⁾,…,𝐸^(𝑁), of dimension 𝐼₁,…,𝐼_𝑁. Tensor 𝒜 has rank one if, for all 𝑛=1,…,𝑁, 𝑁 vectors 𝐮^(𝑛)∈𝐸^(𝑛) exist, such as 𝒜 is equal to their external product: 𝒜=𝐮⁽¹⁾∘⋯∘𝐮^(𝑁). (11)

By extension, 𝒜 is a tensor of rank 𝐾, if 𝐾 is the minimal number or tensor of rank 1 which, by summation, give 𝒜. We write the rank of tensor 𝒜 by Rank(𝒜)=𝐾.

𝑛-Mode Rank
It is possible to define 𝑛-mode rank of a tensor as the generalization of column vectors lines vectors rank of a matrix. 𝑛-mode rank of a tensor 𝒜∈ℂ^{𝐼1×⋯×𝐼𝑁}, written Rank_𝑛(𝒜), is the rank of flattened matrix 𝐀_𝑛 of 𝒜 in 𝑛-mode, Rank_𝑛(𝒜)=Rank(𝐀_𝑛).

Rank-(𝐾₁,…,𝐾_𝑁)
A tensor 𝒜∈ℂ^{𝐼1×⋯×𝐼𝑁} has rank-(𝐾₁,…,𝐾_𝑁) if, for all 𝑛=1,…,𝑁, Rank_𝑛(𝒜)=𝐾_𝑛.

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 𝑁=4 which corresponds to the fourth-order tensor of the TCC values, 𝒯𝒞𝒞∈ℂ^{𝐼1×𝐼2×𝐼3×𝐼4}. The corresponding 𝑛-mode rank values are denoted 𝐾₁,𝐾₂,𝐾₃,𝐾₄. We assume that the value of 𝐾_𝑛 is known or is estimated (see Section 4.2 to know how to estimate them) for all 𝑛=1 to 4.

As an extension to the vector and matrix cases, in the tensor formulation, the projectors on the 𝑛th-mode vector spaces are determined by computing the rank-(𝐾₁,𝐾₂,𝐾₃,𝐾₄) approximation of 𝒯𝒞𝒞 in the least-squares sense. From a mathematical point of view, the rank-(𝐾₁,𝐾₂,𝐾₃,𝐾₄) approximation of 𝒯𝒞𝒞 is represented by tensor 𝒯𝒞𝒞 which minimizes the quadratic tensor Frobenius norm ‖𝒯𝒞𝒞−𝒯𝒞𝒞‖² subject to the constraint that 𝒯𝒞𝒞∈ℂ^{𝐼1×𝐼2×𝐼3×𝐼4} is a rank-(𝐾₁,𝐾₂,𝐾₃,𝐾₄) tensor. The TUCKER3 tensor decomposition also known as Higher-Order SVD (HOSVD) [15] and lower rank-(𝐾₁,𝐾₂,𝐾₃,𝐾₄) tensor approximation (LRTA-(𝐾₁,𝐾₂,𝐾₃,𝐾₄)) 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) The least square solution implies that𝒟=𝒯𝒞𝒞×₁𝐔^(1)𝑇×₂𝐔^(2)𝑇×₃𝐔^(3)𝑇×₄𝐔^(4)𝑇. (13)

Thus the best lower approximation of tensor 𝒯𝒞𝒞 denoted 𝒯𝒞𝒞 is to orthogonally project, for every 𝑛-mode, the vectors of tensor 𝒯𝒞𝒞 on the 𝑛th-mode useful subspace 𝐸₁^(𝑛), for all 𝑛=1 to 4 [9, 13]𝒯𝒞𝒞=𝒯𝒞𝒞×₁𝐏⁽¹⁾×₂𝐏⁽²⁾×₃𝐏⁽³⁾×₄𝐏⁽⁴⁾. (14)

In (14), 𝐏^(𝑛)=𝐔^(𝑛)𝐔^(𝑛)𝑇 is called the projector of the 𝑛th-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 𝑛-mode useful subspaces. The description of TUCKALS3 algorithm, used in rank-(𝐾₁,𝐾₂,𝐾₃,𝐾₄) approximation, is provided in Algorithm 1.

(1) Input: data tensor 𝒯𝒞𝒞, and dimensions 𝐾1,𝐾2,𝐾3,𝐾4 of all 𝑛-mode signal subspaces. (2) Initialization  𝑘=0: For 𝑛=1 to 4, calculate the projectors 𝐏0(𝑛) given by HOSVD-(𝐾1,𝐾2,𝐾3,𝐾4):  (a) 𝑛-mode unfold 𝒯𝒞𝒞 into matrix 𝐓𝐂𝐂𝑛;  (b) Compute the SVD of 𝐓𝐂𝐂𝑛;  (c) Compute matrix 𝐔0(𝑛) formed by the 𝐾𝑛 eigenvectors associated with the 𝐾𝑛 largest singular values of     𝐓𝐂𝐂𝑛. 𝐔0(𝑛) is the initial matrix of the 𝑛-mode signal subspace orthogonal basis vectors;  (d) Form the initial orthogonal projector 𝐏0(𝑛)=𝐔0(𝑛)𝐔(𝑛)𝑇0 on the 𝑛-mode signal subspace;  (e) Compute the HOSVD-(𝐾1,𝐾2,𝐾3,𝐾4) of tensor 𝒯𝒞𝒞 given by   ℬ0=𝒯𝒞𝒞×1𝐏0(1)×2⋯×4𝐏0(4); (3) ALS loop: Repeat until convergence, that is, for example, while ‖ℬ𝑘+1−ℬ𝑘‖2>𝜀, 𝜀>0 being a prior    fixed threshold,  (a) For 𝑛=1 to 4:   (i) Form ℬ(𝑛),𝑘: ℬ(𝑛),𝑘=𝒯𝒞𝒞×1𝐏(1)𝑘+1×2𝐏(2)𝑘+1×3𝐏(3)𝑘+1×4𝐏(4)𝑘+1;   (ii) 𝑛-mode unfold tensor ℬ(𝑛),𝑘 into matrix 𝐁𝑛(𝑛),𝑘;   (iii) Compute matrix 𝐂(𝑛),𝑘=𝐁𝑛(𝑛),𝑘𝐓𝐂𝐂𝑇𝑛;   (iv) Compute matrix 𝐔(𝑛)𝑘+1 composed of the 𝐾𝑛 eigenvectors associated with the 𝐾𝑛 largest eigenvalues      of 𝐂(𝑛),𝑘. 𝐔𝑘(𝑛) is the matrix of the 𝑛-mode signal subspace orthogonal basis vectors at the 𝑘th      iteration;   (v) Compute 𝐏(𝑛)𝑘+1=𝐔(𝑛)𝑘+1𝐔(𝑛)𝑇𝑘+1;  (b) Compute ℬ𝑘+1=𝒯𝒞𝒞×1𝐏(1)𝑘+1×2⋯×4𝐏(4)𝑘+1;  (c) Increment 𝑘. (4) Output: the estimated signal tensor is obtained through 𝒯𝒞𝒞=𝒯𝒞𝒞×1𝐏𝑘(1)stop×2⋯×4𝐏𝑘(4)stop. 𝒯𝒞𝒞 is    the rank-(𝐾1,𝐾2,𝐾3,𝐾4) approximation of 𝒯𝒞𝒞, where 𝑘stop is the index of the last iteration after the    convergence of TUCKALS3 algorithm.

Algorithm 1 

Lower rank-(𝐾₁,𝐾₂,𝐾₃,𝐾₄) tensor approximation (LRTA-(𝐾₁,𝐾₂,𝐾₃,𝐾₄))

In this algorithm, the second-order statistics comes from the SVD of matrix 𝐓𝐂𝐂_𝑛 at step 2b, which is equivalent, up to 1/𝑀_𝑛 multiplicative factor, to the estimation of tensor 𝒯𝒞𝒞𝑛-mode vectors [9]. The definition of 𝑀_𝑛 is given in (8). In the same way, at step 3(a)iii, matrix 𝐂^(𝑛),𝑘 is, up to 1/𝑀_𝑛 multiplicative factor, the estimation of the covariance matrix between tensor 𝒯𝒞𝒞 and tensor ℬ^(𝑛),𝑘𝑛-mode vectors. According to step 3(a)i, ℬ^(𝑛),𝑘 represents data tensor 𝒯𝒞𝒞 filtered in every 𝑚th-mode but the 𝑛th-mode, by projection-filters 𝐏_𝑙^(𝑚), with 𝑚≠𝑛, 𝑙=𝑘 if 𝑚>𝑛 and 𝑙=𝑘+1 if 𝑚<𝑛. 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-(𝐾₁,𝐾₂,𝐾₃,𝐾₄) approximation can simply be achieved by computing the HOSVD-(𝐾₁,𝐾₂,𝐾₃,𝐾₄) of tensor 𝒯𝒞𝒞 [9, 15]. Indeed, the HOSVD-(𝐾₁,𝐾₂,𝐾₃,𝐾₄) of 𝒯𝒞𝒞 consists of the initialization step of TUCKALS3 algorithm, and hence can be considered as a suboptimal solution for the rank-(𝐾₁,𝐾₂,𝐾₃,𝐾₄) 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 the 𝐾 orthonormal basis vectors is to use Gram-Schmidt method (Algorithm 2).

(1) Choose 𝐾, the number of principal axes or eigenvectors required to estimate. Consider matrix 𝐓 and set    𝑝←1. (2) Initialize eigenvector 𝐮𝑝 of size 𝑑×1, for example, randomly; (3) Update 𝐮𝑝 as 𝐮𝑝←𝐓𝐮𝑝; (4) Do the Gram-Schmidt orthogonalization process 𝐮𝑝←𝐮𝑝−∑𝑗=𝑝−1𝑗=1(𝐮𝑇𝑝𝐮𝑗)𝐮𝑗; (5) Normalize 𝐮𝑝 by dividing it by its norm: 𝐮𝑝←𝐮𝑝/‖𝐮𝑝‖. (6) If 𝐮𝑝 has not converged, go back to step 3. (7) Increment counter 𝑝←𝑝+1 and go to step 2 until 𝑝 equals 𝐾.

Algorithm 2 

Fixed-point

The eigenvector with dominant eigenvalue will be measured first. Similarly, all the remaining 𝐾−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 (𝑝−1)th basis vectors will be used to find the 𝑝^th basis vector. The algorithm for 𝑝th basis vector will converge when the new value 𝐮⁺_𝑝 and old value 𝐮_𝑝 are such that 𝐮_𝑝^+𝑇𝐮_𝑝=1. It is usually economical to use a finite tolerance error to satisfy the convergence criterion ‖𝐮_𝑝^+𝑇𝐮_𝑝−1‖<𝜂, where 𝜂 is a prior fixed threshold.

Let 𝐔=[𝐮₁,𝐮₂,…,𝐮_𝐾] be the matrix whose columns are the 𝐾 orthonormal basis vectors. Then 𝐔𝐔^𝑇 is the projector onto the subspace spanned by the 𝐾 eigenvectors associated with the K largest eigenvalues. This subspace is also called “signal subspace.” It can be used in LRTA-(𝐾₁,𝐾₂,𝐾₃,𝐾₄) to retrieve the basis vectors 𝐔₀^(𝑛) in step 2c of Algorithm 1. Thus, the initialization step is faster since it does not need the 𝐼_𝑛 basis vectors but only the 𝐾_𝑛 first ones and it does not need the step 2b, that is SVD of the data tensor 𝑛-mode unfolding matrix 𝐓𝐂𝐂_𝑛.

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 each 𝑛-mode rank. A method to determine them is developed in next subsection.

4.2. 𝑛-Mode Rank Estimation

Whereas a common way to obtain 𝐾₁,𝐾₂,𝐾₃,𝐾₄ values is to use empirical data [18], we wish to present here a scientific tool optimizing the calculation. Each projector 𝐏^(𝑛), 𝑛=1,2,3,4, is estimated by the truncation of unfolding matrices 𝐓𝐂𝐂_𝑛 using SVD, that is, by keeping the 𝐾_𝑛 eigenvectors associated with the 𝐾_𝑛 largest singular values of 𝐓𝐂𝐂_𝑛,𝑛=1,2,3,4. In order to estimate the value of 𝐾_𝑛 for each 𝑛-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 𝑛-mode unfolding matrix 𝐓𝐂𝐂_𝑛. Thus, this method presents the advantage of being self-consistent as it computes the optimal number of kernels for any illumination configuration.

Consequently, for each 𝑛-mode unfolding of 𝒯𝒞𝒞, the form of detection criterion AIC can be expressed asAIC(𝑘)=−2(𝐼_𝑛−𝑘)𝑀_𝑛ln𝒢(𝜆_𝑘+1,…,𝜆_𝐼𝑛)𝒜(𝜆_𝑘+1,…,𝜆_𝐼𝑛)+2𝑘(2𝐼_𝑛−𝑘), (15) and the MDL criterion is given byMDL(𝑘)=−(𝐼_𝑛−𝑘)𝑀_𝑛ln𝒢(𝜆_𝑘+1,…,𝜆_𝐼𝑛)𝒜(𝜆_𝑘+1,…,𝜆_𝐼𝑛)+12𝑘(2𝐼_𝑛−𝑘)ln𝑀_𝑛, (16) where 𝒢 and 𝒜 are geometricand arithmetic mean, respectively. (𝜆_𝑖)_{1≤𝑖≤𝐼𝑛} are 𝐼_𝑛 eigenvalues of the data matrix of the 𝑛-mode unfolding 𝒯𝒞𝒞: 𝜆₁≥𝜆₂≥⋯≥𝜆_𝐼𝑛, and 𝑀_𝑛 is the number of columns of the 𝑛-mode unfolding 𝒯𝒞𝒞.

The 𝑛-mode rank 𝐾_𝑛 is the value of 𝑘 (𝑘=1,…,𝐼_𝑛−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. Figure 3 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 Figure 3.

Figure 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

Figure 4(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.

(a) AIC and MDL estimation criteria

(b) Normalized energy of matrix eigenvalues

Figure 4 

Energy of matrix eigenvalues and comparison between AIC and MDL criteria.

5.3. Fixed Point Algorithm Reconstruction Error

In order to assess our algorithm precision toward SVD decomposition, Figure 5 shows reconstruction error of both fixed point and SVD algorithms, ‖𝐓−𝐓^′‖/‖𝐓‖, where 𝐓 denotes the original matrix and 𝐓^′ 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

Figure 6 shows a comparison of SOCS and LRTA algorithms reconstruction error and both algorithms runtime.

(a) Comparison of LRTA and SOCS reconstruction error

(b) Comparison between SOCS and LRTA runtime

Figure 6 

Comparison between LRTA and SOCS.

To 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 Section 5.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 with Matlab on a 2.4 GHz dual core Pentium with 4Go RAM under Windows XP.

5.5. Comparison of Aerial Images

Figure 7 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.

(a)

(b)

(c)

Figure 7 

Aerial image comparison between SOCS, LRTA, and full Hopkins.

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 Figure 8 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.

Acknowledgment

The authors would like to thank the anonymous reviewers for their careful reading and their fruitful remarks, which have contributed in improving the quality of the paper.