Open Access
Add to BibSonomy
Abstract
We implement a limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) optimization of phase-only computer-generated hologram for a multi-depth three-dimensional (3D) target. Instead of computing the full 3D reconstruction of the hologram, we use a novel method using L-BFGS with sequential slicing (SS) for partial evaluation of the hologram during optimization that only computes loss for a single slice of the reconstruction at every iteration. We demonstrate that its ability to record curvature information enables L-BFGS to have good quality imbalance suppression under the SS technique.
Published by Optica Publishing Group under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
1. INTRODUCTION
Holography, taking its name from the Greek word $o\lambda o\sigma$ (holos), meaning whole, was first introduced in 1948 by Dennis Gabor [1] and was originally named as wavefront reconstruction [2]. It is a technology that can record and reconstruct the wavefront of three-dimensional (3D) objects. Similar to two-dimensional (2D) photography, the earliest holography uses a piece of film to record the diffraction pattern of an object, which can then reconstruct the wavefront showing that object. In order to generate a hologram for objects that do not physically exist, computer-generated holography (CGH) emerged, where a hologram can be calculated through various algorithmic approaches and then be displayed on a spatial light modulator (SLM) modulating a coherent light source in order to reconstruct target images [3–6], either 2D or 3D, which can provide full depth cues at arbitrary angles instead of stereoscopic displays which need to recompute the left-eye and right-eye images every time the position changes. So the multi-depth reconstruction ability is a major advantage of holography technology.
Currently available SLMs can only modulate either phase or amplitude, so algorithms are needed to compute amplitude-only or phase-only holograms, among which the phase-only holograms are usually preferred due to their higher energy efficiency. The classic algorithms include direct binary search [4], simulated annealing [7], and Gerchberg–Saxton [6], which is still widely used despite its being nearly 50 years old. With the developments in modern numerical optimization methods and computation power, advances in CGH algorithms can be made. A literature review has found some recent work that compute CGH using numerical optimization methods [8–12], but speed and quality are still the major challenges in multi-depth hologram generation. The most common multi-depth CGH optimization methods either evaluate their loss against the entire multi-depth 3D target, which is time-consuming, or evaluate the hologram for each plane and then sum the holograms, which introduces quality degradation.
This paper therefore extends on the previous research, which proved the ability of the limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithm to generate a phase-only hologram for a 2D image [13] and proposes the sequential slicing (SS) technique for the optimization of CGH for a multi-depth 3D target, which evaluates the loss for a single slice at each iteration, aiming for quicker hologram generation with proper overall quality and low-quality imbalance across the multiple depths enabled by the second-order nature of the L-BFGS optimization algorithm. The following sections start from the background knowledge of numerical optimization including L-BFGS algorithm, then introduces and carries out the optimization of multi-depth CGH with the SS technique, with results analyzed.
2. METHOD
A. L-BFGS Optimization Algorithm Background
Numerical optimization is a very useful tool to find an optimal solution ${{\bf x}^*}$ that minimizes an objective function $f({\bf x})$. Any unconstrained optimization problem can be written as
and the iteration is given by Eq. (2), where the positive scalar ${\alpha _k}$ is the step length and the vector ${{\bf p}_k}$ is the search direction [14]. There are several optimization algorithms to determine the search direction ${{\bf p}_k}$. The widely used gradient descent (GD) algorithm, as its name suggests, directly uses the negative of the gradient as the search direction [15],In case the first derivative cannot be evaluated directly, it can be approximated using Eq. (4) [14],
However, GD is a greedy first-order optimization method that only considers the current iteration without any global consideration, so it can be extremely slow on complicated problems [14]. Therefore, Newton’s method emerged as a second-order optimization method, which has a search direction as shown in Eq. (5) [16],
Methods that use the Newton direction have a fast rate of local convergence, typically quadratic [14]. However, in many applications, including the optimization of CGH, the Hessian matrix ${\nabla ^2}f({{\bf x}_k})$ in Eq. (5) cannot be evaluated, and even if it could be evaluated, its inverse would be computationally heavy to calculate. Therefore, a quasi-Newton method of BFGS emerged [16], which approximates the inverse of Hessian iteratively, denoted as ${\boldsymbol H}$, following the process described from Eqs. (6)–(9) [14].
Then, aiming to save memory usage for large-scale problems, the L-BFGS algorithm was invented [17], which is listed in Algorithm 1. The L-BFGS algorithm stores the Hessian matrix implicitly and computes the search direction from the $m$ most recent iterations. The L-BFGS algorithm is a quasi-Newtonian method that determines the gradient with curvature information [17] from the gradient history.
Algorithm 1. L-BFGS Algorithm [17]
The optimization algorithms used for CGH in this paper are GD and L-BFGS. The phase-only constraint of CGH can be easily applied by fixing a constant amplitude of the hologram while keeping its phase varying and being the argument of optimization (${\bf x}$).
B. Hologram Optimization for Multi-Depth Targets
As shown in Fig. 1, the multi-depth target is set up as a collection of $n$ slices (${{\bf T}_1}$ to ${{\bf T}_n}$); each slice ${{\bf T}_i}$ is at a distance ${z_i}$ to the hologram plane. And the reconstruction for the hologram ${\bf H}$ at each distance ${z_i}$ is computed using the Fresnel diffraction formula, as shown in Eq. (10),
Fig. 1. Loss between the multi-depth targets (${{\bf T}_1}$ to ${{\bf T}_n}$) and the reconstructions (${{\bf R}_1}$ to ${{\bf R}_n}$) of hologram ${\bf H}$.
Then, to create an objective function [$f({\bf x})$ in Eq. (1)] to optimize, we need to quantify the difference between each target slice (${{\bf T}_i}$) and the respective reconstruction (${{\bf R}_i}$) numerically, and then minimize such difference. Among the many options available, the loss functions selected are mean squared error (MSE) [18], cross entropy (CE) [19], and relative entropy (RE) [20]. To adapt the loss functions for the 2D target image ${{\bf T}_i}$ and reconstructed image ${{\bf R}_i}$ of dimension $X \times Y$, the loss functions are adapted as shown in Eqs. (11)–(13),
MSE is adapted as shown in Eq. (11). MSE is a traditional metric averaging the squared error between the target and observed values. CE is adapted as shown in Eq. (12). CE is often used in classification problems, such as language modeling [21]. RE, also called Kullback–Leibler divergence (usually denoted as ${D_{\textit{KL}}}(P||Q)$), is adapted as shown in Eq. (13). For the uniformity of symbols in this paper, relative entropy is denoted as RE. It is a measure of how much a probability distribution $P$ is different from another probability distribution $Q$. Both CE and RE are usually computed between the true probabilistic distribution and the predicted probabilistic distribution. While the images are not probability distributions, the pixel values are normalized to decimal numbers in the range of 0 to 1 so that CE and RE can be applied.
The effectiveness of L-BFGS optimization of phase-only CGH for a single slice target image ($n = 1$) has been demonstrated in the previous research [13]. However, for a 3D target consisted of multiple slices at different depths, the optimization of CGH becomes challenging.
The typical technique is to sum the losses computed for each slice for each iteration during optimization, which is called the sum-of-loss (SoL) method in this paper. At every iteration, it computes the full 3D reconstructions (${{\bf R}_1}, \ldots ,{{\bf R}_n}$) of the hologram ${\bf H}$ at every distance ${z_i}$, requiring a total of $n$ Fourier transforms, and then computes the total loss between all the targe, ${ t}$ slices and the reconstructed slices, therefore fully evaluating the hologram at each step. Although the total number of iterations does not scale with $n$, the number of Fourier transforms performed at each iteration scales up with $n$, as it needs to compute the full multi-depth reconstructions at every iteration
Fig. 3. Layout of the four-slice target (${z_1} = 1\;\text{cm}$, ${z_2} = 2\;\text{cm}$, ${z_3} = 3\;\text{cm}$, ${z_4} = 4\;\text{cm}$).
Fig. 5. Comparison among SS techniques for the four-slice target using (a) GD algorithm, (b) L-BFGS algorithm, (c) GS algorithm, (d) DCGS algorithm; (e) average NMSE; (f) difference between the maximum and minimum NMSE across all slices.
Another universal technique for 3D CGH is to sum the subholograms ${{\bf H}_i}$ generated for each target slice ${{\bf T}_i}$ to form a total hologram based on the principle of superposition, which is called the sum-of-hologram (SoH) method in this paper. For a fixed number of iterations for each subhologram, the total computation scales up linearly with the number of slices $n$. SoH method’s advantage is its ease of implementation: it can compute multi-depth 3D CGH based on any existing single-slice CGH algorithm. Its major disadvantage for phase-only hologram generation is that the final summing of subholograms will result in a nonuniform amplitude hologram, taking the phase of which will result in discarding the amplitude information of the summed hologram, leading to worse reconstructions quality. And also, SoH method suffers from the defocusing effect from one slice to another, causing additional noise,
This paper proposes what we believe is a novel CGH optimization with the SS technique, as shown in the flow chart in Fig. 2, which only computes the loss for a single slice at each iteration (between a reconstruction ${{\bf R}_i}$ at a single distance ${z_i}$ and its according target slice ${{\bf T}_i}$), where $i$ sweeps through the multi-layer 3D target sequentially when the algorithm iterates. When the final slice is reached ($i = n$), it goes back to the first slice ($i \leftarrow 1$). The proposed method only needs to carry out one Fourier transform at each iteration, and the number of iterations does not need to scale up with $n$. So it is expected to be much quicker than SoL and SoH techniques while producing a proper resulting hologram.
3. RESULTS
A. CGH for Four-Slice Targets
The first example 3D target used is consisted of four slices made from the letters “A”, “B”, “C”, and “D”; each has $512 \times 512\;\text{pixels}$. The positions of the four slices range from 1 cm to 4 cm with a 1 cm gap between each other (i.e., ${z_i} = i\;\text{cm}$). The overall layout is shown in Fig. 3.
As there are two optimization algorithms (GD and L-BFGS), three techniques (SoL, SoH, and SS) and three loss functions (MSE, CE, and RE) under consideration, they can form a total of 18 combinations. In order to control the number of variables, all 18 combinations are set to start from the same initial random hologram and run for the same amount of 100 iterations for the optimization of each hologram, on the same laptop (model ASUS ROG Zephyrus M16, which has a CPU of model i7-11800 H and a GPU of model RTX3060). For the L-BFGS algorithm, the gradient history of size 10 ($m = 10$ in Algorithm 1) is used for all techniques and loss functions. And to ensure a sensible comparison, although three different loss functions are used for optimization of CGH, the metric used to assess the quality of multi-depth reconstructions of the final hologram is the normalized mean squared error (NMSE). As there are a total of four slices in this example, the final NMSE of each is computed separately, and the total optimization run time is recorded. The final results are gathered in Fig. 4. As each slice has a different final error, their mean and standard deviation (SD) are computed for investigations. Three columns are color-coded, where green indicates a better result, while red indicates a worse result.
Comparing the mean of final NMSE and the run time of the proposed SS technique to those of the SoL and SoH techniques in Fig. 4, it can be concluded that, for all combinations of optimizers and loss functions attempted, the proposed SS technique runs much more quickly than the existing SoL and SoH techniques, while still providing a proper result, sitting between the SoL and SoH techniques. So the SS technique has not demonstrated an absolute advantage to the SoL technique yet on this four-slice example. Therefore, further investigation is done on a 30-slice example in Section 3.B. However, the SS technique is both quicker and has better reconstruction quality than the SoH technique, demonstrating an absolute advantage. Meanwhile, the combinations with CE as loss function are much slower and have not demonstrated an advantage in the average NMSE; demonstrating an absolute disadvantage, the results with CE loss or SoH method will not be shown in the per-iteration plots.
For comparison among the SS techniques, the NMSE for each slice and their average and maximum difference values are plotted for GD and L-BFGS algorithms with MSE and RE loss functions in Fig. 5, and sequential GS and dynamic compensatory Gerchberg–Saxton (DCGS) [22] are also implemented and plotted for reference purpose. Apart from the L-BFGS algorithm, all the other algorithms are showing a staircase-like plot, where a decrease in error on one slice results in an increase in error on all other slices, and so the final NMSE of each slice is distinguished a lot from the other. The sequential GS algorithm suffers the most from the quality imbalance between each slice, and the sequential GD algorithm follows. The DCGS algorithm benefits from its modification of the inclusion of a weighting factor consisting of historical amplitude; therefore it managed to converge. From the average NMSE plot [Fig. 5(e)], the proposed sequential L-BFGS method does not appear to have the lowest average NMSE, but it has the lowest quality imbalance across the slices, as shown in the maximum difference plot [Fig. 5(f)]. The L-BFGS algorithm mainly benefits from its inclusion of curvature information during optimization, so that the update of hologram ${\bf H}$ at each iteration takes into account not only the loss for the current slice, but also all historical iterations up to the set history size ($m$ in Algorithm 1). So, the NMSE of each slice stays close for L-BFGS, and at each iteration, the NMSE of all slices behaves in the same way, ensuring each slice has a similar quality.
The final holograms and their reconstructions for L-BFGS algorithm with SoL, SoH, and SS techniques are shown in Fig. 6, with GS with SS technique and DCGS also shown as reference. The reconstructed images confirm the SS technique having a quality between the SoH and SoL methods (for the same amount of iterations), and the SS technique has a much better quality imbalance than sequential GS, which has a very clear reconstruction at the fourth slice (letter “D”) because the iteration stopped at the fourth slice but with much worse reconstructions at other slices. Admittedly, the proposed L-BFGS with SS method cannot surpass the GS-based DCGS algorithm yet.
To prove that the proposed method also works for nonbinary targets, another example of a four-slice 3D target was attempted, as shown in Fig. 7, where two of the slices are replaced by an image of a mandrill [23] and an image of a city scene [8], respectively.
As shown in the final holograms and reconstructions in Fig. 8, the proposed L-BFGS with SS technique still managed to converge, with final reconstruction quality sitting between the SoL and SoH methods, and also having a good quality balance across all slices.
B. CGH for a 30-Slice Target
Another example of a 30-slice target sliced from a 3D teapot is attempted for speed analysis when the number of slices goes higher. As shown in Fig. 9, the Utah teapot [24] is sliced into 30 planes, each of 720p resolution ($1280 \times 720$). Each combination in Fig. 4 was run on a laptop (model ASUS ROG Zephyrus M16 with a CPU of model i7-11800 H and a GPU of model RTX3060). The average and maximum difference of NMSE across all slices plotted against time (in Fig. 10 and Fig. 11, respectively) for all combinations (except those with the SoH technique or with CE as loss function, as they have an absolute disadvantage), for clearer plots.
Figure 10 shows that, for both GD and L-BFGS optimization algorithms, the SS technique is faster than the SoL technique. Among the SS techniques, GD, GS, and DCGS algorithms achieved less average NMSE than the proposed L-BFGS with the SS method, but from the maximum difference of the NMSE plot in Fig. 11, it can be shown that the L-BFGS algorithm has less quality imbalance than the GD algorithm and the GS algorithm, although admittedly it is not as good as the DCGS algorithm. Nevertheless, the proposed method has achieved an improvement of phase retrieval in speed and quality imbalance suppression using optimization algorithms.
Fig. 11. Difference between the maximum and minimum NMSE across all slices versus time plot for the 30-slice target.
4. CONCLUSION
This paper has proposed a method using an L-BFGS optimization algorithm to generate a phase-only hologram for a multi-depth target and discussed its suitability with the SS technique. The L-BFGS with the SS method has demonstrated a good suppression on the quality imbalance across the multi-depth slices, benefiting from the nature of being a second-order optimization algorithm, which implicitly records the historical gradients by other slices for the determination of the descent direction. For both GD and L-BFGS optimization algorithms, the SS technique ran faster and produced better reconstruction quality than the simple SoH technique, and it is quicker than the SoH technique when the number of depths gets large. The proposed method has demonstrated great potential for time-constrained optimization of multi-depth CGH.
Funding
Engineering and Physical Sciences Research Council.
Acknowledgment
This work was supported by the Engineering and Physical Sciences Research Council (EPSRC).
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are available in Ref. [25].
REFERENCES
1. D. Gabor, “A new microscopic principle,” Nature 161, 777–778 (1948). [CrossRef]
2. E. Hecht, Optics, 5th ed. (Pearson, 2017).
3. A. J. Cable, “Real-time high-quality two and three-dimensional holographic video projection using the one-step phase retrieval (OSPR) approach,” Ph.D. thesis (University of Cambridge, 2006).
4. M. A. Seldowitz, J. P. Allebach, and D. W. Sweeney, “Synthesis of digital holograms by direct binary search,” Appl. Opt. 26, 2788–2798 (1987). [CrossRef]
5. H. J. Yang, J. S. Cho, and Y. H. Won, “Reduction of reconstruction errors in kinoform CGHS by modified simulated annealing algorithm,” J. Opt. Soc. Korea 13, 92–97 (2009). [CrossRef]
6. R. W. Gerchberg, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik 35, 237–246 (1972).
7. S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by simulated annealing,” Science 220, 671–680 (1983). [CrossRef]
8. J. Zhang, N. Pégard, J. Zhong, H. Adesnik, and L. Waller, “3D computer-generated holography by non-convex optimization,” Optica 4, 1306–1313 (2017). [CrossRef]
9. S. Liu and Y. Takaki, “Optimization of phase-only computer-generated holograms based on the gradient descent method,” Appl. Sci. 10, 4283 (2020). [CrossRef]
10. S. Choi, J. Kim, Y. Peng, and G. Wetzstein, “Optimizing image quality for holographic near-eye displays with Michelson holography,” Optica 8, 143–146 (2021). [CrossRef]
11. C. Chen, B. Lee, N.-N. Li, M. Chae, D. Wang, Q.-H. Wang, and B. Lee, “Multi-depth hologram generation using stochastic gradient descent algorithm with complex loss function,” Optics Express 29, 15089–15103 (2021). [CrossRef]
12. A. Kadis, B. Wetherfield, J. Sha, F. Yang, Y. Wang, and T. D. Wilkinson, “Effect of bit-depth in stochastic gradient descent performance for phase-only computer-generated holography displays,” London Imaging Meeting 3, 36–40 (2022) [CrossRef] .
13. J. Sha, A. Kadis, F. Yang, and T. D. Wilkinson, “Limited-memory BFGS optimisation of phase-only computer-generated hologram for Fraunhofer diffraction,” in Digital Holography and 3-D Imaging, D. Chu, J. Park, C. Cheng, and P. Ferraro, eds., Technical Digest Series (Optica Publishing Group, 2022), paper W3A.3.
14. J. Nocedal and S. J. Wright, Numerical Optimization (Springer, 2006).
15. C. Lemaréchal, “Cauchy and the gradient method,” in Doc Math Extra (2012), pp. 251–254.
16. R. Fletcher, Practical Methods of Optimization (Wiley, 1987).
17. D. C. Liu and J. Nocedal, “On the limited memory BFGS method for large scale optimization,” Math. Programm. 45, 503–528 (1989). [CrossRef]
18. C. Sammut and G. I. Webb, eds., “Mean squared error,” in Encyclopedia of Machine Learning (Springer, 2010), p. 653.
19. G. Cybenko, D. O’Leary, and J. Rissanen, eds., The Mathematics of Information Coding, Extraction and Distribution, Vol. 107 of The IMA Volumes in Mathematics and Its Applications (Springer, 1998).
20. S. Kullback and R. A. Leibler, “On information and sufficiency,” Ann. Math. Statist. 22, 79–86 (1951). [CrossRef]
21. X. Liu, S. Liu, J. Sha, J. Yu, Z. Xu, X. Chen, and H. Meng, “Limited-memory BFGS optimization of recurrent neural network language models for speech recognition,” in 2018 IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP) (2018), pp. 6114–6118.
22. P. Zhou, Y. Li, S. Liu, and Y. Su, “Dynamic compensatory Gerchberg–Saxton algorithm for multiple-plane reconstruction in holographic displays,” Opt. Express 27, 8958–8967 (2019). [CrossRef]
23. A. Weber, “Sipi image database–misc,” 2022, https://sipi.usc.edu/database/database.php?volume=misc.
24. N. Clark, “Utah teapot (solid). stl” (2015), https://commons.wikimedia.org/wiki/File:Utah_teapot_(solid).stl.
25. J. Sha, A. Kadis, F. Yang, Y. Wang, and T. Wilkinson, “Research data supporting ‘multi-depth phase-only hologram optimization using L-BFGS algorithm with sequential slicing’,” University of Cambridge Repository (2022), https://doi.org/10.17863/CAM.89532.
