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 [36]. 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.

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 by where 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 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 where the superscript “” 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 [1214].

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 , and each element . 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

2.3. -mode product

Consider tensor and a matrix with and (natural numbers excluding ) and for all , -mode product between and is defined as Where operator defines -mode product. It generalizes matrix product to tensor and -mode vector product in a particular -mode.

Each element with index of tensor is given by

2.4. Tensorial Scalar Product

Tensorial scalar product is defined as follows. Suppose two order tensors: and . Scalar product between and is given by

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 is defined as a matrix from , with

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

2.7. External Product

External product “” of several vectors belonging, respectively, to , of dimension defines a tensor of rank 1: in such a way that whatever the indices , the element of tensor is defined by product , where is the 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 , . From II-A, tensor is built on vector space: , of dimension . Tensor has rank one if, for all , vectors exist, such as is equal to their external product:

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 .

-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 , written , is the rank of flattened matrix of in -mode, .

A tensor has rank- if, for all , .

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 which corresponds to the fourth-order tensor of the TCC values, . 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 to .

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 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 The least square solution implies that

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 to [9, 13]

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 of all -mode signal subspaces.
(2) Initialization   : For to , calculate the projectors given by HOSVD-( ):
 (a) -mode unfold into matrix ;
 (b) Compute the SVD of ;
 (c) Compute matrix formed by the eigenvectors associated with the largest singular values of
     . is the initial matrix of the -mode signal subspace orthogonal basis vectors;
 (d) Form the initial orthogonal projector on the -mode signal subspace;
 (e) Compute the HOSVD-( ) of tensor given by
(3) ALS loop: Repeat until convergence, that is, for example, while , being a prior
   fixed threshold,
 (a) For to :
  (i) Form : ;
  (ii) -mode unfold tensor into matrix ;
  (iii) Compute matrix ;
  (iv) Compute matrix composed of the eigenvectors associated with the largest eigenvalues
     of . is the matrix of the -mode signal subspace orthogonal basis vectors at the
  (v) Compute ;
 (b) Compute ;
 (c) Increment .
(4) Output: the estimated signal tensor is obtained through . is
   the rank- approximation of , where is the index of the last iteration after the
   convergence of TUCKALS3 algorithm.

In this algorithm, the second-order statistics comes from the SVD of matrix at step 2b, which is equivalent, up to 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 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 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
(2) Initialize eigenvector of size , for example, randomly;
(3) Update as ;
(4) Do the Gram-Schmidt orthogonalization process ;
(5) Normalize by dividing it by its norm: .
(6) If has not converged, go back to step 3.
(7) Increment counter and go to step 2 until equals .

The eigenvector with dominant eigenvalue will be measured first. Similarly, all the remaining basis vectors (orthonormal to the previously measured basis vectors) will be measured one by one in a reducing order of dominance. The previously measured basis vectors will be used to find the basis vector. The algorithm for basis vector will converge when the new value and old value are such that . It is usually economical to use a finite tolerance error to satisfy the convergence criterion 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 , is estimated by the truncation of unfolding matrices using SVD, that is, by keeping the eigenvectors associated with the largest singular values of . 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 as and the MDL criterion is given by where and are geometricand arithmetic mean, respectively. 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 () 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 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 eigenvectors for matrices and up to around eigenvectors for matrices.

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.

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.

5.4. Comparison between LRTA and SOCS Algorithms

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

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

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 . The model used corresponds to a  nm wavelength illumination,  NA and 0.6 , optical diameter is  nm, and grid size is  nm.

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.

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.


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.