CFD results calibration from sparse sensor observations with a case study for indoor thermal map
Abstract
Current CFD calibration work has mainly focused on the CFD model calibration. However no known work has considered the calibration of the CFD results. In this paper, we take inspiration from the image editing problem to develop a methodology to calibrate CFD simulation results based on sparse sensor observations. We formulate the calibration of CFD results as an optimization problem. The cost function consists of two terms. One term guarantees a good local adjustment of the simulation results based on the sparse sensor observations. The other term transmits the adjustment from local regions around sensing locations to the global domain. The proposed method can enhance the CFD simulation results while preserving the overall original profile. An experiment in an airconditioned room was implemented to verify the effectiveness of the proposed method. In the experiment, four sensor observations were used to calibrate a simulated thermal map with data points. The experimental results show that the proposed method is effective and practical.
CFD results calibration, sparse sensor observation, thermal map estimation, image editing, low rank matrix approximation \IEEEpeerreviewmaketitle
1 Introduction
The estimation of the indoor physical field (e.g., thermal map, airflow pattern, and pressure field) is an important problem in airconditioning system design and indoor architecturaldesign. In recent years, the computational fluid dynamics (CFD) simulation has been widely used in building design to predict the thermal map, pollution dispersion, and ventilation performance [1, 2]. The accuracy of the CFD simulation strongly depends on the setting of boundary conditions and CFD models [3]. In practice, the real boundary conditions of the rooms or buildings are unavailable and the idealized version is used instead. Even though CFD modelling has been extensively researched, the reliability of the results remains a key concern when performing CFD simulations.
Currently, in order to produce reliable simulation results, we need to properly set the boundary conditions and adjust the computational model (numerical and/or physical models) input parameters to amend the agreement between the simulated results and the corresponding experimental data [4, 5, 6]. The current CFD model calibration work has mainly focused on minimizing the mismatch between simulated and experimental results. Tens to hundreds of simulation runs are required for the calibration process of finding acceptable model parameters. However, even a well calibrated CFD model cannot remove the gap between simulated and experimental data, which implies that the results of calibrated CFD models still have room for further improvement.
CFD simulations can map the sampled spatial domain to the airflowfield data. It can provide the global information of airflow fields of interest, which is rather expensive or even impossible to obtain from real experimental observations. On the other hand, sensors (e.g., thermocouples) can only practically provide sparse and isolated observations, but the observations are usually much more accurate than the simulation results. Can we calibrate/correct the simulated physical fields with sparse sensor observations? In this paper, we develop a methodology to solve this problem. To the best of our knowledge, no reported publication has focused on this problem.
The recent work [4, 5, 6] and the references therein mainly focused on the calibration of a CFD model. In this paper, we directly calibrate the results of CFD simulations. The proposed methodology can be applied for general CFD simulations and be a valuable additional step of the standard CFD simulation procedure.
In this work, we propose a methodology to calibrate the physical fields obtained from CFD simulations using sparse sensor observations. Since sensor observations are much more accurate than simulation results, we ignore the difference between sensor observations and the real values and use them as the ground truth. In addition, we assume that the simulation results can provide the rough global profile of the real physical fields, i.e., the profile of the simulated physical field is similar to that of the real physical field. The authors have taken inspiration from the image editing technique [7, 8]. As shown in Fig. 1, an image can be globally edited with sparse control samples. The right image is fused by the original image and the sparse control samples. The result of image editing reserves the content of the original image but with the style of the sparse control samples. For the calibration problem of simulated physical field, we can view the sparse sensor observations as the control samples. The calibration problem can be viewed as a simulated physical field editing problem.
We formulated the simulated physical field calibration as an optimization problem. The objective function is a sum of two parts. One part is responsible to narrow the difference between the sensor observations and the simulation results, while the other part guarantees the similarity between the simulated results and the calibrated results. The optimization problem is solved by a numerical method. The proposed method was tested in an airconditioned room. With the help of four sensor observations, indoor thermal map obtained from CFD simulation was well calibrated.
We summarize the contributions of this work as: 1) we proposed a formulation for the calibration of simulated physical fields from sparse sensor observations; 2) we formulated the proposed calibration problem as an optimization problem; 3) we provided a computationally efficient method to solve the optimization problem, and 4) we verified the effectiveness of the proposed methodology in an airconditioned room.
The rest of this paper is organized as follows. In Section II, we formulate the CFDresult calibration problem. In Section III, we detailed the proposed methodology in which we first formulate the calibration problem as an optimization problem, then present a numerical method to solve it, and finally we provide a computationally efficient approximated solution. In Section IV, we calibrated a thermal map of an airconditioned room to verify the effectiveness of the proposed method. The concluding remarks are given in Section VI.
2 Problem statement
We denote a real physical field (e.g., thermal map, airflow pattern, and pressure field) of interest by , where with being the domain of interest. The physical field obtained from CFD simulation is a numerical solution, which can be simply interpolated to be a continuous one. We denote a continuous physical field interpolated from a CFD solution by , which can be described as
(1) 
where represents the error of the simulated physical field. We assume that , and is a small value compared with the magnitude of .
We denote one of the sensor observations of the physical field by
(2) 
where is the measurement noise of the th observation, is the set of all sensing locations, and is the number of available sensor observations. Without loss of generality, we assume that for all .
In practice, and are known. In what follows, we estimate based on and . With the estimated , we can easily correct the simulated physical field .
3 The proposed methodology
3.1 Formulation of the optimization problem
Since for all , we ignore the measurement noise. We denote the error of the simulated physical field at the sensing locations by , and from (1) and (2) we can obtain
In the remainder of this paper, we ignore the measurement noise and view as the real error of .
To estimate the error of the simulated physical field , i.e., , we minimize the following cost function
(3) 
where
(4) 
is a priori known weight, which can be directly calculated from the simulated physical field. Note that represents the affinity between and . Both and are positive variances. It is obvious that and that is very small unless the locations and are near and their corresponding physical value and are similar.
Solving the minimization problem (3), we can find an estimation , which we denote by . With , we can calibrate the simulated physical field , and the calibration result is
(5) 
The cost function (3) is the summation of two parts. The first part is responsible for the regions around sensing locations. In the region around , the value of is forced to be similar to , especially when is similar to . Consequently, in the region around , the physical value is calibrated to be similar to the sensor observation , especially for the subregion in which is similar to . In this way, in the region around , preserves the profile of but the exact values are adjusted to be similar to the sensor observation . The first term in (3) can guarantee a good local calibration for the simulated physical field, but it has insignificant influence on the regions far away from for all , in which .
On the other hand, the second term has a global influence, which can be explained from the following two perspectives: 1) the second term can transmit the adjustment from the local regions around the sensing locations to the whole physical field; 2) the second term can guarantee that in any local region, the gradient of is small if in that region the simulated values are similar, which guarantees the profile of is similar to that of the simulated physical field .
The optimization problem is controlled by three parameters: , , and . , which we call the balance factor, can balance the contribution of the two terms in (3). A small leads to good calibration in the local regions around sensing locations, while a large can globally reduce the gradient of the adjustment at the cost of the calibration accuracy in the local regions around the sensing locations.
is the product of two Gaussian functions. and control the rate of decay of the two gaussian functions, respectively. , which we call the magnitude variance, mainly depends on the dynamic range of the simulated physical field. If the local gradient of the physical field is large, a large is proper. Otherwise the first Gaussian function in (4) will be very small which may dilute the influence of the second Gaussian function. On the other hand, if the thermal map has a low dynamic range, a relatively small is proper.
, which we call the distance variance, controls the scope of the local region around one of the sensing locations in which the physical data can be calibrated by minimizing the first term in (3). In other words, in the particular local region, the sensor observation have a significant influence on the physical field calibration.
3.2 Numerical method to solve the optimization problem
It is difficult for us to directly minimize in (3) to find the estimation of the error, i.e., . Numerical methods can make the optimization problem easily solvable.
We denote the discrete space of the domain in the CFD simulation by , where is the number of mesh points. We assume that , which implies that each sensing location is constrained to be one of the mesh points. Hence, we can denote and , where is the mesh index of the th sensing location. Then, we can approximate the cost function in (3) as
(6) 
where and . We introduce a new matrix , i.e., is the entry of located at the th row and th column.
For a more elegant description,we use the following cost function instead of
(7)  
where , is a vector with all 1 entries, , is the th column of , and is a matrix with all 1 entries in th column and all 0 entries in others. Here, and are diagonal matrix and trace operators, respectively. It is clear that
and
(8)  
where , , , and . Here, is the th row of . Since is symmetric, .
The matrix in (8) is a diagonal matrix, and . It is clear that .
To minimize the cost function , we take the derivative of (7) with respective to and set it as zero, i.e.,
(9) 
From (9), we can easily obtain
(10) 
Substituting into (5), we obtain the calibrated physical field
(11) 
3.3 Computational efficient approximated solution
To avoid directly solving the inverse of the matrix , which has very large dimension, we provide an approximated solution of the inverse problem. A similar method has been previously applied for image editing[7].
Considering the definition of , i.e., , we find that each column of represents the affinities between the point ( whose mesh index is the column index) and the whole physical field. The neighbouring points have similar affinities with other points. Hence, the neighbouring columns have similar entries, especially for the entries far away from the diagonal, which implies that the neighbouring columns of have a large correlation coefficient. Therefore, the matrix has many near zero eigenvalues, and it can be approximated by a low rank matrix.
Since the matrix , can also be approximated by a low rank matrix. We rewrite the symmetric matrix as
(12) 
where , , and .
We assume that the number of dominant eigenvalues of is . Using the eigen decomposition and ignoring all the near zero eigenvalues, we can approximate as
where is an orthogonal matrix, is a diagonal matrix with all eigenvalues of in the diagonal, consists of all dominant eigenvalues of in its diagonal, and has exactly rank . We denote by
where and have the same sizes as and in (12), respectively. Then, we can obtain that
(13) 
Comparing (13) and (12), as illustrated in Fig. 2, we can find
(14) 
Substituting (14) into (12) yields
(15) 
where is the first column of . Using this way, we can find a low rank approximation of , which can both resolve the storage and computational cost problems in solving .
Substituting the approximation of in (15) into (8) and using the Woodbury formula [9], we can obtain
(16) 
4 Experimental verification
4.1 Statement of the experiment
We tested the proposed methodology for a 7.28m3.32m2.5m room, which is shown in Fig. 3. In the room, an Active Chilled Beam system (ACBs) was installed for air conditioning, as shown on the ceiling of the top picture in Fig. 3. One air outlet is at the corner of the room. The walls and the door are made by thermal insulation materials, and the windows are double glazing. Hence, in the CFD simulation, we simply set the walls, the door, and the windows as insulated walls.
As shown in Fig. 3, the room contains four heaters (0.76m1.2m0.02m). As shown in the bottom picture of Fig. 3, the left two heaters are 300W, and the other two are 400W. The height above the ground of the top surfaces of the heaters is . The bottom faces of the heaters were made by insulated materials. One sensor was used to measure the velocity of the inlet cooling air, as shown in the leftmost picture of Fig. 4. Another sensor was used to measure the inlet cooling air temperature, which is shown in the middle picture in Fig. 4. The other eight sensors were used to measure the temperature of a horizontal plane 1.2m above the ground. Four sensors were hung and suspended just above centers of the four heaters, respectively. The four sensors are represented by stars (s1s4) in the bottom picture of Fig. 3. One of the four is shown in the rightmost picture of Fig. 4. Three sensors (s5s7) were placed at the middle plane of the room. The sensor s8 was placed 20cm from the wall. All the eight sensors (s1s8) were 1.2m above the ground.
The experiment was carried out on 26th May, 2016. The four heaters and the aircondition system were turned on at 9:30 am. We collected the velocity and temperature of the inlet cooling air (via the sensors and , respectively) from 11:17 to 13:10. The data is shown in Fig. 5. The sensors s1s7 recorded the temperature of the 1.2m plane from 9:30 to 13:10, the data for which is shown in Fig. 6. The sensor s8 does not store the sensor readings. We have recorded that after 11:00, the readings of sensor s8 remained at .
4.2 CFD simulation
The simulation was performed using the commercial CFD software Ansys FLUENT v.13.0. The proposed calibration methodology has the potential to be applied for dynamic CFD simulations if both the realtime simulation results and sensor observations are available. In this work, however, for simplicity we only consider the steady state simulation results.
The indoor air was assumed as the ideal gas with density of . The viscous model was chosen as the realizable model. Since the testbed was built in a laboratory, we did not consider the solar radiation. The door, ground, ceiling, windows, and walls are assumed to be insulated walls. During the experiment, all lights were turned off. For the velocity and temperature of the inlet cooling air, we use the mean values of the sensor readings in the time period 12:10 to 13:10 in Fig. 5. The bottom faces of the heaters were made by insulated material, but they were not completely insulated. We assumed that of energy was transferred via the bottom faces of the heaters. Corresponding to the four sensors s1s4, we denote the four heaters by heater 1, heater 2, heater 3 and heater 4, respectively. The boundary conditions of the CFD simulation are detailed in Table 1.
Boundary  Type  Heat transfer  Mass & momentum 

The two inlets  velocityinlet  
Center outlet  pressureoutlet  
Corner outlet  pressureoutlet  
Walls, windows, ceiling, & ground  wall  Insulated  No slip wall 
Top and side faces of heaters 1&2  wall  No slip wall  
Bottom faces of heaters 1&2  wall  No slip wall  
Top and side faces of heaters 3&4  wall  No slip wall  
Bottom faces of heaters 3&4  wall  No slip wall 

inlet cooling air temperature; inlet cooling air speed; the angle of with the ceiling; outlet air temperature; pressure relative to the reference pressure (); convective heat flux.
We created the structured mesh using Gambit v2.4. Since all the walls are insulated, we did not consider the boundary layer of the walls and uniformly divided the horizontal plane ( plane) into squares. In the direction, we made the mesh dense near the ceiling () and heater top faces (), and divided the entire 2.5m height into 137 intervals. The total mesh size is .
To verify the grid independency, we made the mesh dense and increased the mesh size to . We did CFD simulations for both mesh files. The results were found to be in good agreement. Therefore, the mesh file with 8,342,103 cells is proper and the number of cells is large enough for our experiment.
4.3 The calibration results
We use the sensor observations obtained from sensors s1s4 to calibrate the thermal map in Fig. 10. The four sensors s1s4 were just above the centers of the four heaters, respectively.
We set and . The thermal map in Fig. 10 consists of 60955 temperature values of all the mesh points at the plane m. The balance factor is related to the number of mesh points and the number of sensor observations. We set
(18) 
where . The first part in (6) is the summation of terms while the second part is the summation of terms. If we set , it means that we have the same trust on the sensor observations and the simulation results on all mesh points. Obviously, sensor observations are more accurate. We set and the calibrated thermal map is shown in Fig. 10. Here, implies that the influence of each sensor observation in (6) is 100 times important than the simulated data on each mesh point.
Comparing Fig. 10 with Fig. 10, we can easily find that the calibrated thermal map preserves the profile of the CFD results. However, Fig.10 cannot show how the simulated thermal map is calibrated. We show the estimated error (i.e., ) in Fig. 11, which is the difference between Fig. 10 and Fig. 10. Considering (11), we can obtain the simulated thermal map by subtracting the estimated error from the simulated thermal map.
4.4 The influence of balance factor
To check the effect of the balance factor, we set as 0.01, 0.1, 0.5, and 1 respectively. The corresponding estimated errors of the thermal map in Fig. 10 is given in Fig 11. The observations of sensors s1s4 are used to estimate the error. To show the effectiveness of the proposed method, we use the observations of sensor s5s8 as the ground truth to test the estimation performance. The error of the CFD simulation results and the calibrated thermal map at the 8 sensing locations are shown in Fig. 12.
Fig. 11 shows that with different balance factors, the estimated errors (i.e., ) almost have the same profiles. However, the magnitude of the estimated errors decrease when the balance factor increases, which implies that the influence of sensor observations decreases and we can only provide a very small adjustment for the original simulated thermal map. In other words, the balance factor can control the magnitude of the estimated errors, i.e., the magnitude of the adjustment for the simulated thermal map. As shown in Fig. 11 and Fig. 12, with a smaller balance factor, the magnitude of the estimated error is larger and we can obtain a better performance.
Data at the sensing locations  RMSE(s1s8)  improvement(s1s8)  RMSE(s5s8)  improvement(s5s8) 

CFD simulation result  1.1197  /  1.3518  / 
calibrated result ()  0.8536  23.77%  1.0999  18.55% 
calibrated result ()  0.7996  28.59%  1.0451  22.59% 
calibrated result ()  0.7045  37.08%  0.9301  31.06% 
calibrated result ()  0.6794  39.32%  0.8985  33.38% 
 Note: Both and are fixed. We set and . Here,
.
The root mean square error (RMSE) of the thermal map data at sensing locations are given in Table 2. Both Fig. 12 and Table 2 show that we can enhance the thermal map obtained from CFD simulation with the four sensor observations (s1s4). With the calibration, the RMSE of the data at sensing locations s5s8 can be improved by . However, we can see from Fig. 12 that the errors of the regions around s6 and s7 are still large.
The CFD simulation results around s6 and s7 have large errors, i.e., around and , respectively. The proposed method reduces the errors to around and , which are still large. To address this issue, we should consider the capacity of the proposed method. The largest magnitude to which the proposed method can adjust the CFD results almost equals the maximum error of the CFD results at all sensing locations where the sensor observations were used to calibrate the CFD results. It can be seen in Fig. 12 that the maximum error of the CFD results at sensing locations s1s4 is around . Therefore, the maximum value of the estimated error of the simulated thermal map is less than . The balance factor can control the magnitude of the maximum value of the estimated errors. In other words, the balance factor can control the magnitude of the adjustment for the simulated thermal map, but the upper bound of this magnitude is determined by the sensing locations.
4.5 The influence of magnitude variance
To check the influence of magnitude variance, we set as 1000, 100, 10, and 1, respectively. With the four sensor observations (s1s4), we estimated the error of the thermal map in Fig. 10. The estimated errors are given in Fig.13. The error of the estimated thermal map at the 8 sensing locations are shown in Fig. 14, and the RMSE of the calibrated thermal map at the 8 sensing locations are given in Table 3.
Fig. 13 shows that with different magnitude variances, the maximum magnitude of the estimated errors (i.e., ) are almost the same, which implies that different with the balance factor, the magnitude variance has no significant influence on the the maximum magnitude of the estimated errors.
It can be seen from Fig. 1313 that the region in the two red circle is nearer to sensor s4 than s3 but the estimated error is similar with the estimated error around s3. We can see form Fig. 10 that the temperature values of the thermal map in the region marked by the two circles (in Fig. 1313) are similar to those in the region around sensor s3. The proposed thermal map calibration work is a process of fusing sensor information and the information of CFD results. A smaller magnitude variance can enhance the influence of the CFD results in this information fusion work.
In addition, Fig. 13 shows that if in the regions far away the sensing locations s1s4, the estimated errors are almost zero. In other words, the proposed method has no influence for these regions. As we mentioned before, a small magnitude variance leads to a small value of the first Gaussian function in (4), which can dilute the influence of the second Gaussian function. The smaller the magnitude variance , the smaller the weight , especially when is far away the sensing location . Therefore, if the magnitude variance is too small, the proposed method could not globally calibrated the thermal map.
Data at the sensing locations  RMSE(s1s8)  improvement(s1s8)  RMSE(s5s8)  improvement(s5s8) 

CFD simulation result  1.1197  /  1.3518  / 
calibrated result ()  1.0142  9.42%  1.4339  6.05% 
calibrated result ()  0.7652  31.66%  1.0543  21.91% 
calibrated result ()  0.6957  37.87%  0.9267  31.31% 
calibrated result ()  0.6794  39.32%  0.8985  33.38% 
 Note: Both and are fixed. We set and . Here,
.
Fig. 14 shows that with the decrease of the magnitude variance, the errors around sensing locations s1s4 also decrease but those around sensing locations s5s8 do not. From Table 3 we can find that the RMSE at sensing locations s5s8 increases with respect to the decrease of the magnitude variance. If we reduce the magnitude variance to 1, the calibration work makes the accuracy of the simulated thermal map worse. Hence, we conclude that reducing the magnitude variance can improve the local calibration but has no significant help for global calibration. If enough sensor observations are available we can set a small magnitude variance unless a large one is proper.
4.6 The influence of distance variance
To check the influence of distance variance, we set as 1, 0.5, 0.2, and 0.1, respectively. With the four sensor observations (s1s4), we estimated the error of the thermal map in Fig. 10. The estimated errors are given in Fig.15. The error of the estimated thermal map at the 8 sensing locations are shown in Fig. 16, and the RMSE of the calibrated thermal map at the 8 sensing locations are given in Table 4.
Data at the sensing locations  RMSE(s1s8)  improvement(s1s8)  RMSE(s5s8)  improvement(s5s8) 

CFD simulation result  1.1197  /  1.3518  / 
calibrated result ()  0.9558  14.64%  1.3518  0% 
calibrated result ()  0.9594  14.32%  1.3568  0.37% 
calibrated result ()  0.7117  36.44%  1.0062  25.45% 
calibrated result ()  0.6794  39.32%  0.8985  33.38% 
 Note: Both and are fixed. We set and . Here,
.
From Fig. 15 we can find that like the magnitude variance, the distance variance has no significant influence on the maximum value of the estimated errors. It is clear shown in Fig. 15 that with a smaller distance variance, the scope of nonzero region of the estimated errors dramatically decrease. Fig. 16 shows that when the distance variance , the proposed method can provide almost perfect calibration at the four sensing locations (s1s4). Table 4 shows that with a smaller distance variance, the errors of the calibrated thermal map at the sensing locations s5s8 increase. These agree with our previous analysis that the distance variance controls the scope of the region on which the sensor observation has significant influence. A small distance variance can lead to very good local calibration while a large can provide global influence for an isolated sensor observation.
5 concluding remarks
CFD simulation results calibration from sparse sensor observations is a new and interesting problem. We formulated this problem as an optimization problem. The proposed method can provide an effective solution for this problem. The experimental results show that with four sensor observations (s1s4), the proposed method has improvement on the accuracy of the simulation data around sensors s5s8. If we consider the whole thermal map (including the region around sensors s1s4), the improvement can be even larger.
We need to set three parameters for the proposed method: balance factor, magnitude variance, and distance variance. The balance factor controls the magnitude of the adjustment for the simulated thermal map, but the upper bound of this magnitude is determined by the sensing locations. The distance variance controls the scope of the regions on which an isolated sensor observation has significant influence. A small distance variance can lead to very good local calibration while a large one can provide global influence for an isolated sensor observation. The magnitude variance controls the influence of the CFD results and scales the influence of the distance variance. Reducing the magnitude variance can improve the local calibration but has no significant help for global calibration. If enough sensor observations are available we can set both small magnitude variance and small distance variance. If the number of available sensor observations is very limited, we should set relatively larger values for the two variances to achieve global performance.
In addition, the capacity and the performance of the proposed method is closely related to the number of sensor observations and the sensing locations. The issue of how many sensors are required and where to place them are very interesting topics for future research.
Acknowledgements
This work was supported by Singapore’s National Research Foundation under NRFCRP8201103 and partially supported by the Energy Research Institute at NTU (ERI@N). The first author would like to thank Mr. Ke Ji, Mr. Xingyu Zhang, and Mr. Huynh Nam Khoa for their help in the experiments.
References
 Z. Zhai, “Application of computational fluid dynamics in building design: aspects and trends,” Indoor and Built Environment, vol. 15, no. 4, pp. 305–313, 2006.
 Q. Chen, “Ventilation performance prediction for buildings: A method overview and recent applications,” Building and Environment, vol. 44, no. 4, pp. 848–858, 2009.
 J. Srebric, V. Vukovic, G. He, and X. Yang, “CFD boundary conditions for contaminant dispersion, heat transfer and airflow simulations around human occupants in indoor environments,” Building and Environment, vol. 43, no. 3, pp. 294–303, 2008.
 M. Hajdukiewicz, M. Geron, and M. M. Keane, “Formal calibration methodology for CFD models of naturally ventilated indoor environments,” Building and Environment, vol. 59, pp. 290–302, 2013.
 S. Guillas, N. Glover, and L. MalkiEpshtein, “Bayesian calibration of the constants of the k– turbulence model for a CFD model of street canyon flow,” Computer Methods in Applied Mechanics and Engineering, vol. 279, pp. 536–553, 2014.
 O. T. Kajero, R. B. Thorpe, T. Chen, B. Wang, and Y. Yao, “Kriging metamodel assisted calibration of computational fluid dynamics models,” AIChE Journal, 2016.
 X. An and F. Pellacini, “Appprop: allpairs appearancespace edit propagation,” in ACM Transactions on Graphics (TOG), vol. 27, no. 3. ACM, 2008, p. 40.
 L. Xu, Q. Yan, and J. Jia, “A sparse control model for image and video editing,” ACM Transactions on Graphics (TOG), vol. 32, no. 6, p. 197, 2013.
 G. H. Golub and C. F. Van Loan, Matrix computations. JHU Press, 2013.