Abstract
In this paper, we propose to adapt the multilinear algebra tools to the tensor of Transmission CrossCoefficients (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 wellknown RETs is the optical proximity correction (OPC), which consists of making mask precompensation 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. ModelBased 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 fourway 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 CrossCoefficients (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 fourway array. will denote the size of tensor in dimension . Tensor approach permits to treat a fourway 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 twodimensional 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 socalled โuseful subspace.โ Useful subspacebased methods are applied to source localization in array processing and image denoising. These methods were adopted to fourwayarray (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 Nway 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 zeroorder tensor is a scalar, a firstorder tensor is a vector, and a secondorder 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 thmode vector space of dimension , associated with the thmode of tensor . By definition, is generated by the column vectors of the thmode flattening matrix. The thmode 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 thmode vectors of tensor . The lower rank tensor approximation is possible by extending the classical subspace approach to tensors by assuming that, whatever the thmode 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 wellknown (Akaike Information Criterion) AIC or (Minimum Description Length) MDL criteria [10]. These are entropybased 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 IIA, 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, .
Rank
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 fourthorder 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 thmode vector spaces are determined by computing the rank() approximation of in the leastsquares 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 HigherOrder 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 thmode useful subspace , for all to [9, 13]
In (14), is called the projector of the thmode 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.

In this algorithm, the secondorder 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 thmode but the thmode, by projectionfilters , 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 HOSVDbased 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 GramSchmidt 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 GramSchmidt method (Algorithm 2).

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 wellknown 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 selfconsistent 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 fourway 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 fourthorder 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
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.
(a) Comparison of LRTA and SOCS reconstruction error
(b) Comparison between SOCS and LRTA 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 Tshaped 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)
(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 Figure 8 using the same calculus only with 4 kernels for Lshaped 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 CrossCoefficients 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 fourway 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.