We study a real-time bidding problem resulting from a set of
contractual obligations stipulating that a firm win a specified
number of heterogenous impressions or ad placements over a defined duration in a realtime
auction.  The contracts specify item targeting criteria (which may
be overlapping), and a supply requirement.
Using the Pontryagin maximum principle, we show that the resulting
continuous time and time inhomogenous planning problem can be
reduced into a finite dimensional convex optimization problem and
solved to optimality.  In addition, we provide algorithms to update
the bidding plan over time via a receding horizon.  Finally, we
provide numerical results based on real data and show a connection
to production-transportation problems.

Computational Advertising; Realtime Bidding; Optimal
Control; Auction Theory; Second Price Auction; Production
Transportation Problem

We acknowledge the support of the Natural Sciences and
Engineering Research Council of Canada (NSERC), [funding reference
number 518418-2018].  Cette recherche a été financée par le Conseil
de recherches en sciences naturelles et en génie du Canada (CRSNG),
[numéro de référence 518418-2018].

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

Eindhoven University of Technology Eindhoven

titleQueueing Systems

journaltitleManagement Science
titleRepeated auctions with budgets in ad exchanges: Approximations and design

Taylor & Francis

journaltitleJournal of the Operational Research Society
titleNonlinear programming

booktitleProceedings of the Tenth ACM International Conference on Web Search and Data Mining
titleReal-time bidding by reinforcement learning in display advertising

New York, NY, USA

titleModel Predictive Control—Classical, Robust and Stochastic

booktitleProceedings of the Eighth International Workshop on Data Mining for Online Advertising
titleA dynamic pricing model for unifying programmatic guarantee and real-time bidding in display advertising

journaltitleInformation Systems Research
titleOnline display advertising markets: A literature review and future directions

Springer Science & Business Media

titleFunctional analysis, calculus of variations and optimal control

booktitleProceedings of the 17th ACM SIGKDD international conference on Knowledge discovery and data mining
titleBid landscape forecasting in online ad exchange marketplace

journaltitleExpositiones Mathematicae
titleOn the substitution rule for Lebesgue–Stieltjes integrals

journaltitlearXiv preprint arXiv:2001.06587
titleScalable Bid Landscape Forecasting in Real-time Bidding

booktitleProceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining
titleRecurrent Neural Networks for Stochastic Control in Real-Time Bidding

journaltitleAvailable at SSRN 2066175
titleOptimal bidding strategies and equilibria in dynamic auctions with budget constraints

titleIAB internet advertising revenue report conducted by PricewaterhouseCoopers (PWC)

ACM New York, NY, USA

journaltitleACM SIGMETRICS Performance Evaluation Review
titleBidding with limited statistical knowledge in online auctions

booktitle2014 American Control Conference
titleAdaptive control using Heisenberg bidding

journaltitleEuropean Journal of Control
titleControl problems in online advertising and benefits of randomized bidding strategies

journaltitleTransportation Science
titleThe transportation-production problem

booktitleProceedings of the Eighth International Workshop on Data Mining for Online Advertising
titleiPinYou global rtb bidding algorithm competition dataset

journaltitlearXiv preprint arXiv:1204.0535
titleDoubleclick ad exchange auction

Morgan & Claypool Publishers

journaltitleSynthesis Lectures on Communication Networks
titlePerformance modeling, loss networks, and statistical multiplexing

journaltitleEconometrica: Journal of the Econometric Society
titleA decomposition algorithm for solving the multifacility production-transportation problem with nonlinear production costs

New York, NY, USA

Association for Computing Machinery

journaltitleProc. ACM Meas. Anal. Comput. Syst.
titleAn Optimal Bidding Algorithm for Online Ad Auctions with Overlapping Targeting Criteria

journaltitleOnline: http://cvxopt. org/documentation/coneprog. pdf
titleThe CVXOPT linear and quadratic cone program solvers

Wiley Online Library

journaltitleThe Journal of finance
titleCounterspeculation, auctions, and competitive sealed tenders
family=van der Walt,
family=Jarrod Millman,
family=Vand erPlas,
family=van Mulbregt,

journaltitleNature Methods
titleSciPy 1.0: Fundamental Algorithms for Scientific Computing in Python

Now Publishers, Inc.

journaltitleFoundations and Trends® in Information Retrieval
titleDisplay advertising with real-time bidding (RTB) and behavioural targeting

booktitleJoint European Conference on Machine Learning and Knowledge Discovery in Databases
titleFunctional bid landscape forecasting for display advertising

Springer Science & Business Media

titleAll of nonparametric statistics

journaltitlearXiv preprint arXiv:1209.4340
titleMoments and absolute moments of the normal distribution

booktitleProceedings of the 27th ACM International Conference on Information and Knowledge Management
titleBudget constrained bidding by model-free reinforcement learning in display advertising

booktitleProceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining
titlePredicting winning price in real time bidding with censored data

booktitleProceedings of the Seventh International Workshop on Data Mining for Online Advertising
titleReal-time bidding for online advertising: measurement and analysis

journaltitleManagement Science
titleSoft floors in auctions

booktitleProceedings of the Ninth ACM International Conference on Web Search and Data Mining
titleFeedback control of real-time display advertising

booktitleProceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining
titleOptimal real-time bidding for display advertising

journaltitlearXiv preprint arXiv:1407.7073
titleReal-time bidding benchmarking with ipinyou dataset

booktitleProceedings of the 18th ACM SIGKDD international conference on Knowledge discovery and data mining
titleJoint optimization of bid and budget allocation in sponsored search

booktitleProceedings of the 22nd ACM SIGKDD international conference on Knowledge discovery and data mining
titleBid-aware gradient descent for unbiased learning with censored data in display advertising


# Autoencoder based Randomized Learning of Feedforward Neural Networks for
Regression
††thanks: Supported by Grant 2017/27/B/ST6/01804 from the National Science
Centre, Poland (research) and and Grant 020/RID/2018/19 "Regional Initiative
of Excellence" from the Minister of Science and Higher Education, Poland
(conference).

Grzegorz Dudek  Department of Electrical Engineering
Czestochowa University of Technology
Częstochowa, Poland
<EMAIL_ADDRESS>
###### Abstract

Feedforward neural networks are widely used as universal predictive models to
fit data distribution. Common gradient-based learning, however, suffers from
many drawbacks making the training process ineffective and time-consuming.
Alternative randomized learning does not use gradients but selects hidden node
parameters randomly. This makes the training process extremely fast. However,
the problem in randomized learning is how to determine the random parameters.
A recently proposed method uses autoencoders for unsupervised parameter
learning. This method showed superior performance on classification tasks. In
this work, we apply this method to regression problems, and, finding that it
has some drawbacks, we show how to improve it. We propose a learning method of
autoencoders that controls the produced random weights. We also propose how to
determine the biases of hidden nodes. We empirically compare autoencoder based
learning with other randomized learning methods proposed recently for
regression and find that despite the proposed improvement of the autoencoder
based learning, it does not outperform its competitors in fitting accuracy.
Moreover, the method is much more complex than its competitors.

###### Index Terms:

autoencoder, feedforward neural networks, randomized learning algorithms

##  I Introduction

Feedforward neural networks (FNNs) have attracted a great deal of interest due
to their excellence predictive performance and universal approximation
capabilities. The most popular learning methods involve some kind of gradient
descent algorithm to learn FNN weights iteratively. However, gradient-based
methods suffer from many drawbacks making the learning process ineffective and
time-consuming. This is because they are sensitive to the initial values of
the parameters. The learning trajectory, starting with different initial
parameters, leads to local minima of the loss function. Thus, globally optimal
parameters are not guaranteed. Moreover, the learning process is time-
consuming for complex target functions (TFs), big data sets and large FNN
architectures. Randomized learning, such as a random vector functional link
(RVFL) network [1], has been proposed as an alternative to conventional FNNs
iterative learning using gradients. In randomized learning, the parameters of
the hidden nodes are selected randomly and stay fixed. They do not need to be
tuned during the training stage. The only parameters that need to be learned
are the output weights. This makes the optimization problem convex [2]. As
such, it can be solved easily and quickly using a standard least-squares
method. Despite this simplification, randomized FNN learning still possesses
universal approximation capabilities, provided there are a sufficient number
of nonlinear hidden nodes [1], [3].

Many simulation studies reported in the literature show the high performance
of the randomized FNN when compared to fully adaptable FNNs. Randomization,
which is cheaper than optimization, ensures simplicity of implementation and
faster training. However, a challenging and still open question in randomized
learning is how to choose appropriate weights and biases for the hidden nodes
to ensure best model performance [4], [5]. To deal with this problem, various
methods of generating hidden node parameters have been proposed. The simplest
and most popular solution is to select both weights and biases from a uniform
distribution over some symmetric interval $U=[-u,u]$. Usually, this interval
is assigned as fixed, typically $[-1,1]$, regardless of the data, TF, and type
of activation functions (AFs). The independence of the hidden nodes from data
is seen as an asset. This overly-simplistic approach was criticized as
illogical and misleading [6]. Therefore, to improve its performance,
optimization of interval $U$ for a specified application is recommended. For
example, in [7], a supervisory mechanism which randomly assigns hidden node
parameters from an adaptively selected interval, was proposed. This paper
clearly reveals that the selection of random parameters should be data
dependent to ensure the universal approximation property of the resulting
randomized FNN.

Data dependent random parameters were also recommended in [8]. The authors of
this work noticed that if the hidden nodes are chosen at random and not
subsequently trained, they are usually not placed in accordance with the
density of the input data. In such a case, training of linear parameters is
less ineffective at reducing errors. Therefore, in order to improve learning
performance, the authors advise unsupervised placement of hidden nodes
according to the input data density. This recommendation was implemented in
the methods proposed in [9] and [10]. In these works, it was noticed that as
the weights and biases of hidden nodes have different functions, they should
not be selected from the same interval. The weights decide about AF slopes and
should reflect TF complexity, while the biases decide about the placement of
AF in the input space. The biases should ensure the introduction of the most
nonlinear fragments of AFs into the input hypercube. These fragments are most
useful for modeling TF fluctuations. According to the methods described in [9]
and [10], we first select the proper interval for the weights based on the AF
features and TF properties. Then, the biases are calculated based on the
weights and data distribution. In [11], a data-driven method was proposed to
improve further the FNN randomized learning. This method introduces the AFs
into randomly selected regions of the input space and adjusts the slopes of
individual AFs to the TF slopes in these regions. As a result, the AFs mimic
the TF locally and their linear combination approximates smoothly the entire
TF.

An interesting method of generating parameters for RVFL was proposed recently
in [12]. For this, the authors employ support-vector machines which in a
supervised manner, by solving their corresponding optimization problems,
generate the pre-trained weights. These weights are used to initialize the
hidden layer of the proposed RVFL architecture. An alternative approach to
generating hidden nodes in FNNs is unsupervised parameter learning using
autoencoders (AEs), which was first introduced in [13]. AE, in the encoding
phase, transforms input data into a meaningful feature representation obtained
from the hidden layer. Then, in the decoding phase, this feature
representation is converted to the original inputs. The information hidden in
original data can be explored and encoded into the output weights of AE [14].
These output weights are then introduced to FNN as hidden node weights instead
of randomly generated weights [13]. This approach was applied in [15] for
classification tasks. Here, RVFL uses a sparse AE with $\ell_{1}$-norm
regularization to adaptively learn superior hidden node parameters for
specific learning tasks. The authors claim that the learned network parameters
in their sparse pre-trained RVFL are embedded with the valuable information
about input data, which alleviates the randomly generated parameter issue and
improves algorithmic performance.

Another classifier based on RVFL with unsupervised parameter learning was
proposed in [16]. In this solution, randomization based stacked AEs with a
denoising criterion are used to extract better, higher-level representations.
Each randomization based AE acts as an independent feature extractor and a
deep network is obtained by stacking several such AEs. The network is built
hierarchically with high level feature extraction followed by a final
classification layer, which is RVFL with direct links.

The authors of both works on AE based randomized learning, [15] and [16],
report experimental results on many real-world classification data sets from
different domains. The results confirm the excellent effectiveness of the
proposed solutions. Encouraged by these state-of-the-art results for FNN
classifiers trained in an unsupervised manner using AEs, in this study, we
analyze unsupervised parameter learning of FNN for regression. We compare this
approach with alternative methods, which were proposed recently in [10].

This paper makes the following contributions:

  1. 1.

We analyze AE based unsupervised parameter learning for a FNN regression model
and find that this method has some drawbacks. We show how to improve it. We
propose a learning method of AE that controls the produced random weights for
FNN. We also propose how to determine the biases for FNN.

  2. 2.

We empirically compare AE based learning with other randomized learning
methods proposed recently for regression and find that AE based learning does
not outperform its competitors in fitting accuracy, and, in fact, it is much
more complex.

The remainder of this paper is structured as follows. In Section II, we
describe randomized FNN learning and methods for generating random parameters.
AE based generating of random parameters is presented and critically analyzed
from the perspective of AF distribution and shaping in Section III. In Section
IV, we analyze the complexity of randomized AE. The performance of AE based
randomized learning is evaluated in Section V. Finally, in Section VI, we
conclude the work.

##  II Randomized Learning of FNN

Let us consider a shallow FNN architecture with $n$ inputs, a single-hidden
layer including $m$ nonlinear nodes, and a single output. AFs of hidden nodes,
$h_{i}(\mathbf{x})$, map nonlinearly input vectors
$\mathbf{x}=[x_{1},x_{2},...,x_{n}]^{T}$ into $m$-dimensional feature space.
An output node combines linearly $m$ nonlinear transformations of the inputs.
FNN expresses a function in the form:

$\varphi(\mathbf{x})=\sum_{i=1}^{m}\beta_{i}h_{i}(\mathbf{x})$ (1)

where $\beta_{i}$ is the output weight between the $i$-th hidden node and the
output node.

Such shallow architecture has a universal approximation property, even when
the hidden layer parameters are not trained but generated randomly from the
proper distribution [1], [3].

The output weights
$\boldsymbol{\beta}=[\beta_{1},\beta_{2},...,\beta_{m}]^{T}$ can be determined
by solving the following linear problem:
$\mathbf{H}\boldsymbol{\beta}=\mathbf{Y}$, where
$\mathbf{H}=[\mathbf{h}(\mathbf{x}_{1}),\mathbf{h}(\mathbf{x}_{2}),...,\mathbf{h}(\mathbf{x}_{N})]^{T}\in\mathbb{R}^{N\times
m}$ is the hidden layer output matrix, and
$\mathbf{Y}=[y_{1},y_{2},...,y_{N}]^{T}$ is a vector of target outputs. Using
Moore-Penrose pseudoinverse, the optimal solution is given by:

$\boldsymbol{\beta}=\mathbf{H}^{+}\mathbf{Y}$ (2)

where $\mathbf{H}^{+}$ denotes the Moore–Penrose generalized inverse of matrix
$\mathbf{H}$.

The hidden node parameters, i.e. weights
$\mathbf{a}=[a_{1},a_{2},...,a_{n}]^{T}$ and bias $b$, control AF slope,
orientation and position in the input space. For a sigmoid AF given by the
formula:

$h(\mathbf{x})=\frac{1}{1+\exp\left(-\left(\mathbf{a}^{T}\mathbf{x}+b\right)\right)}$
(3)

weight $a_{j}$ expresses the sigmoid slope in the $j$-th direction and bias
$b$ decides about the sigmoid shift along a hyperplane containing all
$x$-axes. An appropriate selection of the slopes and shifts of all sigmoids
determine the approximation properties of the model.

As it was shown in [9] and [10], the standard way of selecting both the hidden
node weights and biases randomly, from the same interval, $a_{i,j},b_{i}\sim
U(-u,u)$, is misguided. This is because the optimal interval for weights,
which determine the AF slope range, is not the optimal interval for biases,
which represent an AF shift. And vice versa. With this in mind, in [9] and
[10] separate methods for determining weights and biases were proposed. The
approach described in [9] first selects the weights $a_{i,j}$ from $U(-u,u)$.
The bounds of the interval, $u$, are adjusted to TF complexity. For flat TFs
we expect lower bounds, while for strongly fluctuating TFs we expect higher
bounds. Once the weights have been selected, the biases are determined in such
a way as to ensure the steepest fragments of the sigmoids (which are around
their inflection points) are introduced into input hypercube
$H=[x_{1,\min},x_{1,\max}]\times...\times[x_{n,\min},x_{n,\max}]$. The
resulting equation for the $i$-th hidden node bias is as follows:

$b_{i}=-\mathbf{a}_{i}^{T}\mathbf{x}_{i}^{*}$ (4)

where $\mathbf{x}_{i}^{*}=[x_{i,1}^{*},...,x_{i,n}^{*}]$ is a point from $H$
where the $i$-th sigmoid has one of its inflection points.

As you can see from (4), the biases are dependent on the weights. Point
$\mathbf{x}_{i}^{*}$ can be selected as follows:

  * •

this can be some point randomly selected from $H$: $\mathbf{x}_{i}^{*}\sim
U(H)$. This method is suitable when the input points are evenly distributed in
$H$.

  * •

this can be some randomly selected training point:
$\mathbf{x}_{i}^{*}=\mathbf{x}_{\xi}\in\Phi$, where $\xi\sim U\\{1,...,N\\}$,
and $\Phi$ is a training set. This method distributes the sigmoids according
to the data density, avoiding empty regions.

  * •

this can be a prototype of the training point cluster:
$\mathbf{x}_{i}^{*}=\mathbf{p}_{i}$, where $\mathbf{p}_{i}$ is a prototype of
the $i$-th cluster. This method requires the clustering of training points
into $m=$#nodes clusters.

It was noticed in [10] that the relationship between weights $a$ and the slope
angles of sigmoids $\alpha$ is highly nonlinear. The standard interval for
$a$, $U=[-1,1]$ corresponds to the interval
$U_{\alpha}=[-14^{\circ},14^{\circ}]$ for $\alpha$, so only flat sigmoids are
obtainable in such a case. To get steep sigmoids, with $\alpha$ near
$90^{\circ}$, the bounds for $U$ should be $u>100$. For narrow $U$, such as
$[-1,1]$, the distribution of $\alpha$ is similar to a uniform one. When the
interval for $a$ is extended, the $\alpha$ distribution changes such that
larger angles, near the bounds of $U_{\alpha}$, are more probable than smaller
ones. When $a\in[-100,100]$, more than $77\%$ of sigmoids are inclined at an
angle greater than $80^{\circ}$, so they are very steep. In such a case, there
is a real threat of overfitting. To generate sigmoids with uniformly
distributed slope angles, first, we select randomly $|\alpha_{i,j}|\sim
U(\alpha_{\min},\alpha_{\max})$, where the bound angles,
$\alpha_{\min}\in(0^{\circ},90^{\circ})$ and
$\alpha_{\max}\in(\alpha_{\min},90^{\circ})$, are adjusted to the TF
complexity. Then, we calculate the weights from:

$a_{i,j}=4\tan\alpha_{i,j}$ (5)

Finally, to introduce the sigmoids into the input hypercube $H$, the biases
are calculated from (4).

Both methods of generating random parameters of hidden nodes described above,
we use in the experimental part of the work as comparative methods for AE
based method. We denote them as R$a$M, i.e. a random weights $a$ method, and
R$\alpha$M, i.e. a random slope angles $\alpha$ method.

Fig. 1 illustrates an approximation of a highly nonlinear TF by FNN trained
using R$\alpha$M (similar results were obtained for R$a$M). TF, shown as the
dashed line in the left panel, is fitted accurately by the function built by
FNN (red line). This fitted function is composed of 25 hidden node sigmoids,
which are shown in the right panel. Note that the steepest fragments of the
sigmoids are introduced by R$\alpha$M into the input interval, which is shown
by the gray field in the panel on the right. This fragments are the most
useful for modeling the TF fluctuations. The saturated AF fragments in the
input interval are avoided. The interval for sigmoid slope angles,
$U_{\alpha}=[0^{\circ},83^{\circ}]$, was adjusted to the TF complexity. This
gives a perfect fitting.

Figure 1: Fitted curve (left panel) and hidden node sigmoids (right panel) for
R$\alpha$M with $|\alpha|\sim U(0^{\circ},83^{\circ})$.

##  III Autoencoder based Generation of Hidden Nodes Parameters

Unsupervised parameter learning using AEs was proposed in [13]. In this work,
randomization based autoencoders (RAE) are used for unsupervised feature
extraction for a multilayer FNN classifier. RAE consists of two parts, an
encoder and a decoder. The encoder maps the input data randomly into some
latent representation in such a way that the decoder is able to reconstruct
the original input data form this representation. RAE is a single hidden layer
FNN with $n$ inputs, $n$ outputs, and $m$ nonlinear hidden nodes. The sigmoid
AFs, $g(\mathbf{x})$, are used for the hidden layer and linear AFs are used
for the output layer. The sigmoid AF is given by:

$g(\mathbf{x})=\frac{1}{1+\exp\left(-\left(\mathbf{w}^{T}\mathbf{x}+c\right)\right)}$
(6)

where $\mathbf{w}=[w_{1},w_{2},...,w_{n}]^{T}$ are the hidden node weights and
$c$ is its bias.

In the first step, the learning method using RAE (RAEM) selects randomly the
hidden node parameters for RAE, i.e. weights $\mathbf{w}_{i}$ and biases
$c_{i},i=1,2,...,m$. Typically, both are taken from the uniform distribution
and interval $[-1,1]$. In the second step, the hidden layer output matrix is
calculated:
$\mathbf{G}=[\mathbf{g}(\mathbf{x}_{1}),\mathbf{g}(\mathbf{x}_{2}),...,\mathbf{g}(\mathbf{x}_{N})]^{T}\in\mathbb{R}^{N\times
m}$, where
$\mathbf{g}(\mathbf{x}_{l})=[g_{1}(\mathbf{x}_{l}),g_{2}(\mathbf{x}_{l}),...,g_{m}(\mathbf{x}_{l})]$
is a vector of hidden node outputs for input pattern $\mathbf{x}_{l}$.
Finally, the output weight matrix,
$\mathbf{V}=[\mathbf{v}_{1},\mathbf{v}_{2},...,\mathbf{v}_{m}]^{T}\in\mathbb{R}^{m\times
n}$, where $\mathbf{v}_{i}=[v_{i,1},v_{i,2},...,v_{i,n}]$, is calculated from:

$\mathbf{V}=\mathbf{G}^{+}\mathbf{X}$ (7)

where $\mathbf{G}^{+}$ denotes the Moore–Penrose generalized inverse of matrix
$\mathbf{G}$ and $\mathbf{X}\in\mathbb{R}^{N\times n}$ is an input matrix.

Due to randomized learning of RAE, the output weights, $\mathbf{V}$, can be
obtained easily using a standard linear least-squares method. These weights
are considered to be the latent features of input data [13], and are used as
the hidden node weights for FNN instead of random weights. Thus,
$\mathbf{A}=\mathbf{V}^{T}$, i.e. $a_{j,i}=v_{i,j}$.

As for the biases of hidden nodes $b_{i}$, it is hard to find in the
literature, how they are selected in RAE based FNN learning. The exception is
[15], where the authors outline the way of which the biases were determined.
They calculate the bias for the $i$-th node as the mean value of this node
weights: $b_{i}=1/n\sum_{j=1}^{n}a_{j,i}$. Let us look at this case from the
perspective of AF distribution in the input space.

When the bias is defined as the weight average, the inflection point of the
sigmoid for the one-dimensional case, which is for $h(x)=0.5$, can be obtained
from:

$h(x)=\frac{1}{1+\exp\left(-\left(ax+b\right)\right)}=0.5\rightarrow
x=-\frac{b}{a}=-1$ (8)

where we substituted $b=a$ assuming that a bias is the mean value of the
weights.

Thus, the inflection points of all sigmoids are in $x=-1$. This case is
visualized in Fig. 2. Note that the steepest fragments of the sigmoids, which
are around point $x=-1$, are far from the input interval. This interval
includes many saturated fragments of the sigmoids, which results in poor
fitting. The accuracy does not improve with the number of hidden nodes.

Figure 2: Fitted curve (left panel) and hidden node sigmoids (right panel) for
RAEM with $b_{i}=a_{i}$.

For the $n$-dimensional case, the sigmoid value in the inflection points is:

$h(\mathbf{x})=\frac{1}{1+\exp\left(-\left(\mathbf{a}^{T}\mathbf{x}+b\right)\right)}=0.5$
(9)

Substituting the mean value of weights for $b$ in this equation, after
transformations we obtain:

$\mathbf{a}^{T}\mathbf{x}+\overline{a}=0\\\ $ (10)

where $\overline{a}=\frac{1}{n}\sum_{j=1}^{n}a_{j}$.

This equation expresses the inflection hyperplane of the sigmoid. As we can
see, this hyperplane is totally dependent on weights $a_{j}$. So, the weights
generated by the RAE determine all sigmoid features, slopes in all directions
and shift. The interdependence of the slopes and shift is an undoubted
disadvantage of this approach. It is unjustified and limits the model’s
flexibility. The slopes should correspond to the TF complexity and the shift
should be related to data distribution. Unfortunately, RAEM proposed in [15]
cannot control separately the slopes and shifts of the sigmoids when it
generates them. When ignoring biases and assuming $b_{i}=0$ for all hidden
nodes, all sigmoids pass through $\mathbf{x}=(0,...,0)$. This is also an
unacceptable solution as it limits the approximation properties of the model.

In an attempt to improve the RAEM performance, we use the same solution for
biases as in R$a$M and R$\alpha$M. We calculate them from (4) on the basis of
weights $\mathbf{a}_{i}$ produced by RAE and randomly selected points
$\mathbf{x}^{*}$. Thus, we distribute the sigmoids in $H$ according to the
data distribution. The result for the one-dimensional case is shown in Fig. 3.
As we can see, the accuracy did not improve significantly. Moreover, when we
add new hidden nodes (up to 2000), the RMSE remains unacceptable, i.e. above
0.1. This means that the problem is not only due to the mis-determination of
the biases, but also due to the too flat sigmoids that do not correspond to TF
complexity.

Figure 3: Fitted curve (left panel) and hidden node sigmoids (right panel) for
RAEM with $b_{i}=-a_{i}x^{*}_{i}$.

The sigmoid slopes in FNN are determined by the RAE output weights, $v$. These
weights are dependent on the RAE random projection $\mathbf{G}$, which in turn
is dependent on the RAE hidden node random parameters: $w$ and $c$. In the
above described simulations, these parameters were both selected from the
standard interval $U_{AE}=[-1,1]$. As we shown in Section II, this is an
incorrect approach. Thus, to improve the regression model performance, we
propose to use R$a$M for generating $w$ and $c$. According to this method, we
optimize the interval for $w$, $U_{AE}=[-u_{AE},u_{AE}]$, and calculate the
biases $c$ analogously to (4):

$c=-\mathbf{w}^{T}\mathbf{x}^{*}$ (11)

As a result of this modified RAEM learning, we expect that RAE will produce
the appropriate weights $\mathbf{A}=\mathbf{V}^{T}$, which provide FNN
sigmoids with the slopes adjusted to the TF complexity.

Fig. 4 shows the effect of the interval $U_{AE}$ bounds on the median of
absolute values of weights $v$ (left panel) and on the fitting error (right
panel) for TF shown in Figs. 1-3 (the number of hidden nodes was $m=25$,
results are averaged over 100 runs). Note that median of $|v|$ decreases
quickly with $u_{AE}$. RMSE reaches its minimum for $u_{AE}\in(0.05,0.2)$. For
such $U_{AE}$ bounds, RAE produces weights $v$ whose $median|v|$ is in the
range from around 5 to 14 (for comparison, when weights $w$ were selected from
the standard interval $U_{AE}=[-1,1]$, $median|v|$ was around 1.5). For such a
case, RMSE reaches an acceptable level of 0.005, which is similar to those
obtained by R$a$M and R$\alpha$M for the same number of hidden nodes.

Fig. 5 depicts fitting results for $U_{AE}=[-0.1,0.1]$. Note steeper sigmoids
than in the cases presented in Fig. 3, resulting in good fitting. It should be
noted that the optimal interval $U_{AE}$ depends on the number of hidden
nodes. When, instead of 25, we used 200 hidden nodes, the optimal values of
the interval bounds were $u_{AE}\in(0.005,0.022)$. In such a case $median|v|$
was from around 4 to 17. Compare these values with $median(|v|)\approx 0.19$
which we obtain for $U_{AE}=[-1,1]$.

Figure 4: The effect of the interval $U_{AE}$ bounds on the median of absolute
values of weights $v$ (left panel) and on the fitting error (right panel).

Figure 5: Fitted curve (left panel) and hidden node sigmoids (right panel) for
RAEM with $w_{i}\in[-0.1,0.1]$ and $c_{i}=-w_{i}x^{*}_{i}$.

In the above analysis, RAE was trained without regularization. However,
$\ell_{1}$, $\ell_{2}$ and elastic-net regularization are widely used in RAE
to prevent overfitting [13], [15], [16]. Regularization in RAE decreases
weights $v$ to improve the generalization property of the model, i.e.
generalization of mapping $\mathbf{x}$ to themselves. From the point of view
of the regression FNN model trained using RAEM, regularisation in RAE is an
unfavorable operation because it further flattens FNN sigmoids. As we showed
above, it is a disadvantage for strongly fluctuating TFs.

Note that RAEM in its standard version (without optimization of $U_{AE}$)
produces weights $v$ dependent only on the distribution of points $\mathbf{x}$
in the input space and independently of the TF. So, for two TFs with different
complexity (one of them flat and the second one with strong fluctuations),
having the same distributed x-points, RAE can produce exactly the same weights
$v$. This must be considered a serious drawback.

These findings on RAEM can be summarised in tree points:

  * •

the output weights $v$ produced by RAE determine the sigmoid slopes in the FNN
regression model and are crucial for accurate TF approximation. These weights
are dependent on the RAE hidden nodes parameters, i.e. weights $w$ and biases
$c$. The standard way of generating both these parameters from the same
interval $U_{AE}$ is unjustified and misleading. We recommend optimizing this
interval for weights $w$ and determining biases $c$ from (11).

  * •

the optimal interval $U_{AE}$ for $w$ depends on the number of hidden nodes
$m$. Thus, for each value of $m$ considered, interval $U_{AE}$ should be
optimized.

  * •

the hidden nodes biases of the FNN regression model, $b$, determine the
sigmoid placement in the input space. They should introduce the steepest
fragments of the sigmoids into the input hypercube. Incorrectly selected, such
as $b_{i}=\overline{a}_{i}$ or $b_{i}=0$, they lead to the placement of the
saturation parts of the sigmoids into $H$. To avoid this, we recommend
determining biases $b$ from (4).

##  IV Complexity of RAE

To generate weights $v$, RAE can use linear least squares regression, lasso,
ridge regression or elastic net algorithms [16]. In all these cases the total
time complexity is $O(Nm^{2}+m^{3})$ (this is true for lasso when least-angle
regression (LARS) is used for fitting linear regression models [17]). For
$N>m$ the RAE runs in linear time with the size of the training set, $N$, and
in quadratic time with the number of nodes, $m$. For $N\leq m$, it runs in
cubic time with $m$, although in [18] it was shown that for ridge regression
this cost can be reduced to $O(N^{2}m)$.

Taking into account the optimization process, i.e. the selection of
hyperparameters, the time complexity is as follows. In RAE without
regularization the number of hidden nodes $m$ and interval $U_{AE}$ need to be
found using cross-validation. Thus, the computational load increases linearly
with the number of data splits used in cross-validation, $s$, and also with
the number of points of the grid which is $l_{m}l_{u}$, where $l_{m}$ and
$l_{u}$ are the number of searching values for $m$ and $u_{AE}$, respectively.
So, the total complexity will be $O(sl_{m}l_{u}Nm^{2}+sl_{m}l_{u}m^{3})$. In
RAE with $\ell_{1}$ and $\ell_{2}$ regularization, the regularization
parameter $\lambda$ should also be selected in cross-validation. So, the
complexity of RAE with lasso or ridge regression is
$O(sl_{m}l_{u}l_{\lambda}Nm^{2}+sl_{m}l_{u}l_{\lambda}m^{3})$, where
$l_{\lambda}$ is the number of searching points for $\lambda$. Elastic net
combines ridge and lasso, with the tuning parameter $\alpha$ that balances the
weights of ridge against lasso. If the number of searching points for $\alpha$
is $l_{\alpha}$, the total complexity of RAE with elastic net regularization
is
$O(sl_{m}l_{u}l_{\lambda}l_{\alpha}Nm^{2}+sl_{m}l_{u}l_{\lambda}l_{\alpha}m^{3})$.

##  V Simulation Study

In this section, to demonstrate the fitting properties of RAEM, we report some
simulation results over several regression problems. They include an
approximation of extremely nonlinear TFs:

TF1

$g(\mathbf{x})=\sum_{j=1}^{n}\sin\left(20\cdot\exp x_{j}\right)\cdot
x_{j}^{2},\,x_{i}\in[0,1]$

TF2

$g(\mathbf{x})=-{\sum_{i=1}^{n}\sin(x_{i})\sin^{20}\left(\frac{ix_{i}^{2}}{\pi}\right)},\,x_{i}\in[0,\pi]$

TF3

$g(\mathbf{x})=418.9829n-{\sum_{i=1}^{n}x_{i}\sin(\sqrt{|x_{i}|})},\\\
x_{i}\in[-500,500]$

We considered these functions with $n=2,5$ and 10 arguments. The sizes of the
training and test sets depended on the number of arguments. They were 5000 for
$n=2$, 20,000 for $n=5$, and 50,000 for $n=10$. All arguments for TF3-TF5 were
normalized to $[0,1]$, and the function values were normalized to $[-1,1]$.
Two-argument functions TF1-TF3 are shown in Fig. 6. Note that TF1 combines
flat regions with strongly fluctuating regions, TF2 expresses flat regions
with perpendicular grooves, and TF3 fluctuates strongly, showing the greatest
amplitude at the borders.

Figure 6: Target functions TF1-TF3 for $n=2$.

We use a modified version of RAEM. That is, the interval for the hidden nodes
in RAE, $U_{AE}$, was optimized for each number of hidden nodes. This is
because, as we show in Section III, the optimal interval is dependent on the
hidden node number. To introduce the sigmoids of RAE into input hypercube $H$,
biases $c$ were calculated from (11). And analogously, to introduce the
sigmoids of the FNN regression model into $H$, we calculate biases $b$ from
(4). For comparison, we use R$a$M and R$\alpha$M. For each experiment, we run
100 independent training sessions.

In Table I, we present the fitting test errors (RMSE) of each method for each
TF. The optimal bounds of the intervals from which the model parameters ($w$,
$a$ and $\alpha$, respectively) were randomly selected are also shown. In
R$\alpha$M, we set fixed lower bounds, $\alpha_{\min}=0^{\circ}$, while the
upper bounds, $\alpha_{\max}$, were selected for each TF as $90^{\circ}$. For
RAEM and R$a$M, we observe the difference in the optimal interval sizes
between 2-argument and more than 2-argument TFs. R$a$M provides wide
intervals, i.e. steep sigmoids, for 2-argument TFs, and narrow intervals, i.e.
flat sigmoids, for 5- and 10-argument TFs. Similarly RAEM provides steep
sigmoids for 2-argument TFs (narrow $U_{AE}$ produces higher weights $a$), and
flat sigmoids for 5- and 10-argument TFs. This could be explained by the
change in the TF landscape, which flattens with an increasing number of
dimensions. Interestingly, R$\alpha$M is insensitive to this phenomenon,
giving the same broad interval for $\alpha$ regardless of the number of
arguments.

As can be seen from Table I, R$\alpha$M demonstrates the highest fitting
accuracy for all TFs. This was confirmed by a Wilcoxon signed-rank test with
$\alpha=0.05$. Note that results for RAEM and R$a$M are very similar. Figs.
7-9 depict fitting test errors depending on the number of hidden nodes $m$.
Shaded regions are 10th and 90th percentiles, measured over 100 trials. As can
be seen from these figures, the confidence intervals of RAEM and R$a$M overlap
for higher $m$. This means that these two methods generate similar weights
$a$, which provide a similar set of basis functions for FNN. To study this
issue further, we show in Figs. 10-12 the histograms of weights $a$ generated
by RAEM for all TFs. It is evident from these figures that weight $a$
distributions deviate from the uniform distribution, especially for $n=5$ and
$10$. Thus, RAEM, unlike R$a$M, does not generate uniformly distributed
weights $a$. The weight distributions are unimodal, symmetrical, bell-shaped,
and centered at 0. Note that the intervals for $a$ observed in Figs. 10-12 in
most cases correspond to the intervals $U$ selected as optimal by R$a$M (see
hyperparameter $u$ for R$a$M in Table I).

TABLE I: Results for TF1-TF3. | #nodes | RAEM | R$a$M | R$\alpha$M
---|---|---|---|---
| $m$ | RMSE | $u_{AE}$ | RMSE | $u$ | RMSE | $\alpha_{\max}$
TF1, $n=2$ | 800 | $0.0011\pm 0.00050$ | 0.00100 | $0.0012\pm 0.00093$ | 20 | $\mathbf{0.0006\pm 0.00040}$ | 90
TF2, $n=2$ | 3000 | $0.0051\pm 0.00407$ | 0.00004 | $0.0049\pm 0.00338$ | 100 | $\mathbf{0.0006\pm 0.00051}$ | 90
TF3, $n=2$ | 2000 | $0.0058\pm 0.00095$ | 0.00006 | $0.0059\pm 0.00097$ | 100 | $\mathbf{0.0024\pm 0.00062}$ | 90
TF1, $n=5$ | 200 | $0.2210\pm 0.00012$ | 8 | $0.2219\pm 0.00010$ | 0.1 | $\mathbf{0.2137\pm 0.00374}$ | 90
TF2, $n=5$ | 500 | $0.2381\pm 0.00101$ | 0.1 | $0.2392\pm 0.00012$ | 0.5 | $\mathbf{0.1734\pm 0.00912}$ | 90
TF3, $n=5$ | 500 | $0.2381\pm 0.00051$ | 1 | $0.2392\pm 5\cdot 10^{-8}$ | 0.02 | $\mathbf{0.1717\pm 0.00903}$ | 90
TF1, $n=10$ | 100 | $0.2327\pm 0.00006$ | 6 | $0.2327\pm 0.00005$ | 0.01 | $\mathbf{0.2318\pm 0.00072}$ | 90
TF2, $n=10$ | 100 | $0.2579\pm 0.00043$ | 2 | $0.2577\pm 0.00034$ | 0.1 | $\mathbf{0.2551\pm 0.00290}$ | 90
TF3, $n=10$ | 100 | $0.2237\pm 0.00013$ | 3 | $0.2241\pm 0.00006$ | 0.01 | $\mathbf{0.2188\pm 0.00212}$ | 90

Figure 7: RMSE depending on the hidden node numbers for 2-argument TFs.

Figure 8: RMSE depending on the hidden node numbers for 5-argument TFs.

Figure 9: RMSE depending on the hidden node numbers for 10-argument TFs.

Figure 10: Histograms of weights $a$ generated by RAEM for 2-argument TFs.

Figure 11: Histograms of weights $a$ generated by RAEM for 5-argument TFs.

Figure 12: Histograms of weights $a$ generated by RAEM for 10-argument TFs.

It is clear from Table I and Figs. 7-9 that the best fitting was achieved for
R$\alpha$M. Interestingly, for each TF this method generated weights $a$ from
the same distribution, which is shown in Fig. 13.

Figure 13: Histograms of weights $a$ generated by R$\alpha$M for all TFs.

To compare the performance of RAEM, R$a$M and R$\alpha$M on real-world
regression problems, we performed experiments on several data sets from
different domains. Data sets were collected from the KEEL repository,
http://www.keel.es/ (stock, laser, treasury, dee, and machineCPU data sets)
and from Delve repository,
https://www.cs.toronto.edu/$\sim$delve/data/kin/desc.html (kin8nm data set).
The number of samples and arguments in each data set are shown in Table II.
The input and output variables were normalized into $[0,1]$. The data sets
were divided into training sets containing 75% of the samples selected
randomly, and test sets containing the remaining samples. The optimal values
of hyperparameters, i.e. hidden node numbers and sizes of intervals for random
parameters, were selected by 5-fold cross-validation. Results are shown in
Table II. The bold values indicate the lowest errors while the values in
italics indicate the highest errors (Wilcoxon signed-rank test with
$\alpha=0.05$ was used to compare errors). Note that RAEM demonstrate the
highest errors for four data sets. The best method is R$\alpha$M, which
yielded the lowest errors for five out of six data sets.

TABLE II: Results for real-word regression problems. Data set | #samples/ | RAEM | R$a$M | R$\alpha$M
---|---|---|---|---
| #arguments | RMSE | $m$ | $u_{AE}$ | RMSE | $m$ | $u$ | RMSE | $m$ | $\alpha_{\max}$
stock | 950/9 | $\mathit{0.0303\pm 0.0013}$ | 150 | 0.1585 | $0.0296\pm 0.0016$ | 200 | 4 | $\mathbf{0.0204\pm 0.0006}$ | 200 | 30
laser | 993/4 | $\mathit{0.0193\pm 0.0039}$ | 80 | 0.8913 | $\mathit{0.0204\pm 0.0037}$ | 60 | 0.5 | $\mathbf{0.0160\pm 0.0006}$ | 60 | 45
treasury | 1049/15 | $0.0102\pm 0.0012$ | 150 | 0.3162 | $0.0097\pm 0.0009$ | 100 | 1 | $\mathbf{0.0084\pm 0.0002}$ | 150 | 65
dee | 365/6 | $0.0486\pm 0.0113$ | 10 | 8.9125 | $0.0501\pm 0.0125$ | 10 | 1 | $0.0472\pm 0.0077$ | 10 | 40
machineCPU | 209/6 | $\mathit{0.0486\pm 0.0115}$ | 10 | 7.9433 | $0.0398\pm 0.0104$ | 16 | 0.5 | $\mathbf{0.0290\pm 0.0023}$ | 26 | 60
kin8nm | 8192/8 | $\mathit{0.0654\pm 0.0024}$ | 1000 | 0.0631 | $0.0636\pm 0.0013$ | 1000 | 1 | $\mathbf{0.0514\pm 0.0009}$ | 900 | 20

Fig. 14 compares variants of RAEM:
RAEM1 the improved RAEM proposed in this study (optimized interval for weights
$w$, $U_{AE}=[-u_{AE},u_{AE}]$; biases $c$ and $b$ determined from (11) and
(4), respectivelly),
RAEM2 variant with the fixed interval for $w$, $U_{AE}=[-1,1]$; biases $c$ and
$b$ determined from (11) and (4), respectivelly,
RAEM3 variant with random selection of both $w$ and $c$ from the fixed
interval $[-1,1]$; biases $b$ determined from (4),
RAEM4 variant with random selection of $w$, $c$ and $b$ from the fixed
interval $[-1,1]$,
RAEM5 variant with random selection of both $w$ and $c$ from the fixed
interval $[-1,1]$; biases $b$ determined as mean values of weights,
$b_{i}=\overline{a}_{i}$.

Note that the proposed modification of RAE in most cases gave the lowest
errors compared to other RAE variants. The laser and treasury data sets were
insensitive to the method of generating weights and biases. For these data
sets, the errors for all RAE variants were similar.

Figure 14: Comparison of RAEM variants.

##  VI Conclusion

In this work, we showed that RAE based randomized learning of FNN suffers from
several drawbacks. First, RAE produces the random weights for the FNN
predictive model by taking into account just the input data. This approach is
questionable because the input data does not contain information about TF
complexity. TFs with strong fluctuations need higher weights than flat TFs to
be modeled accurately. Unfortunately, RAE for both these cases can generate
similar sets of weights ignoring completely TF complexity. Second, RAE does
not generate the biases for the FNN predictive model. These biases, which
determine the distribution of the activation functions in the input space, are
crucial for the approximation properties of the predictive model.

In this study, we propose improved, unsupervised parameter learning using
RAEs. First, we introduce the possibility of controlling the magnitude of the
random weights produced by RAE. This is realized by appropriately generating
the RAE hidden node parameters. Second, we determine the biases for the FNN
predictive model so that the sigmoids have their steepest fragments introduced
into an input hypercube. These fragments are the most useful for modeling TF
fluctuations. The proposed modifications make the RAE method more flexible,
more data dependent and more dependent on the complexity of the solved
problem.

The experimental part of the work does not provide evidence that the improved
RAEM outperforms in fitting accuracy other new methods of generating random
parameters. Moreover, its complexity is much greater as it requires additional
learning of RAE. Therefore, applying it to regression problems, rather than
simpler and faster methods, may be questionable.

## References

  * [1] B. Igelnik and Y.H. Pao, "Stochastic choice of basis functions in adaptive function approximation and the functional-link net," IEEE Trans. Neural Network, vol. 6, no. 6, pp. 1320–1329, 1995.
  * [2] J. Principe and B.Chen, "Universal approximation with convex optimization: Gimmick or reality?" IEEE Computational Intelligence Magazine, vol. 10, no. 2, pp. 68–77, 2015.
  * [3] D. Husmeier, "Random vector functional link (RVFL) networks," in Neural Networks for Conditional Probability Estimation: Forecasting Beyond Point Predictions, chap. 6, pp. 87–97. Springer-Verlag, London, 1999.
  * [4] W. Cao, X. Wang, Z. Ming, and J. Gao, "A review on neural networks with random weights," Neurocomputing, vol. 275, pp. 278–287, 2018.
  * [5] L. Zhang and P. Suganthan, "A survey of randomized algorithms for training neural networks," Information Sciences, vol. 364–365, pp. 146–155, 2016.
  * [6] M. Li and D. Wang, "Insights into randomized algorithms for neural networks: Practical issues and common pitfalls," Information Sciences, vol. 382–383, pp. 170–178, 2017.
  * [7] D. Wang and M. Li, "Stochastic configuration networks: Fundamentals and algorithms," IEEE Trans. Cybernetics, vol. 47, no. 10, pp. 3466–3479, 2017.
  * [8] A.N. Gorban , I.Y. Tyukin , D.V. Prokhorov , and K.I. Sofeikov, "Approximation with random bases: Pro- et contra," Information Sciences, vol. 364–365, pp. 129–145, 2016.
  * [9] G. Dudek, "Generating random weights and biases in feedforward neural networks with random hidden nodes," Information Sciences, vol. 481, pp. 33–56, 2019.
  * [10] G. Dudek, "Generating random parameters in feedforward neural networks with random hidden nodes: Drawbacks of the standard method and how to improve it," in: Neural Information Processing. ICONIP 2020. Communications in Computer and Information Science, vol. 1333, pp. 598–606, Springer, Cham, 2020.
  * [11] G. Dudek, "Data-driven randomized learning of feedforward neural networks," 2020 International Joint Conference on Neural Networks (IJCNN), Glasgow, United Kingdom, pp. 1–8, 2020.
  * [12] M. Tanveer, M.A. Ganaie and P.N. Suganthan, "Ensemble of classification models with weighted functional link network," Applied Soft Computing, doi: https://doi.org/10.1016/j.asoc.2021.107322, 2021.
  * [13] L.L.C. Kasun, H. Zhou, G. Huang, C. Vong, "Representational learning with extreme learning machine for big data," IEEE Intelligent Systems, vol. 28, pp. 31–34, 2013.
  * [14] G. Hinton and R.R. Salakhutdinov, "Reducing the dimensionality of data with neural networks," Science, vol. 313, pp. 504–507, 2006.
  * [15] Y. Zhang, J. Wu, Z. Cai, B. Du, and P.S. Yu, "An unsupervised parameter learning model for RVFL neural network," Neural Networks, vol. 112, pp. 85–97, 2019.
  * [16] R. Katuwal and P.N. Suganthan, "Stacked autoencoder based deep random vector functional link neural network for classification,” Applied Soft Computing, vol. 85, 105854, 2019.
  * [17] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani, "Least angle regression," Annals of Statistics, vol. 32, no. 2, pp. 407–499, 2004.
  * [18] T. Hastie, R. Tibshirani, and J. Friedman, The Elements of Statistical Learning: Data Mining, Inference, and Prediction, 2nd ed., Springer-Verlag New York, 2009.


# Enhancement of mechanical squeezing via feedback control

Chao Meng<EMAIL_ADDRESS>Niels Bohr Institute, University of Copenhagen,
2100 Copenhagen, Denmark  Warwick P. Bowen<EMAIL_ADDRESS>Australian Research Council Centre of Excellence in Quantum Biotechnology,
School of Mathematics and Physics, University of Queensland, St Lucia,
Queensland 4072, Australia

###### Abstract

We explore the generation of nonclassical mechanical states by combining
continuous position measurement and feedback control. We find that feedback-
induced spring softening can greatly enhance position squeezing. Conversely,
even with a pure position measurement, we find that spring hardening can
enable momentum squeezing. Beyond enhanced squeezing, we show that feedback
also mitigates degradation introduced by background mechanical modes.
Together, this significantly lowers the barrier to measurement-based
preparation of nonclassical mechanical states at room temperature.

††preprint: APS/123-QED

Measurements play a crucial role in the preparation of quantum states and are
utilized across a variety of fields, including quantum computing [1], quantum
sensing [2, 3], and tests of fundamental physics [4]. Continuous measurement
of the position of a mechanical oscillator has been extensively studied and
has established a standard limit to measurement precision [5], impacting
precision optomechanical sensors and gravitational wave detectors [6, 7].

Recently, it has been shown that a position squeezed mechanical state can be
prepared by continuous measurement as long as the measurement rate is faster
than the rate that noise couples to the position of the oscillator [8]. The
noise coupling rate is related to the mechanical frequency. At room
temperature, this necessitates the use of a low frequency mechanical
resonance. Typical high-frequency mechanical oscillators (above a megahertz)
used in optomechanical experiments [9, 10, 11, 12, 13, 14] cannot enter the
squeezing regime even with state-of-the-art technology. Moreover, background
of mechanical resonances can severely degrade state preparation [15]. This
degradation cannot be evaded by increasing the measurement rate as this also
increases the strength of background interactions. Together, this raises the
question of whether it is practical to generate quantum squeezing at room
temperature via continuous measurement. Furthermore, although momentum
squeezing is crucial for enhancing force sensing, it is not clear how to
achieve this via position measurement.

In this Letter, we show that the above-mentioned limitations can be overcome
through feedback control. We find that feedback-induced spring softening,
shifting the mechanical resonance frequency lower, can improve conditional
position squeezing by a factor equal to the fractional change in resonance
frequency. Moreover, by shifting the resonance away from those of background
mechanical modes, we show that this feedback can greatly mitigate degradation
due to background resonances. Consequently, the criteria for achieving quantum
squeezing are substantially relaxed, facilitating the entry of existing high-
frequency optomechanical devices into the non-classical regime at room
temperature. Strikingly, we also establish that spring hardening, shifting the
mechanical frequency higher, permits momentum squeezing through continuous
position measurement, thus granting optomechanical devices—particularly those
operating at low frequencies with large mass—access to momentum below the
zero-point motion. Collectively, these advancements pave the way for room-
temperature quantum state preparation and quantum sensing.

Figure 1:  Schematic diagram of state-preparation. A laser resonantly excites
an optomechanical cavity. The mechanical position is detected through optical
phase quadrature measurement. Feedback is applied to it to control the
mechanical resonance frequency. The Wiener filters ($H_{x}$ and $H_{p}$)
filter the measurement, providing optimal estimates of $x$ and $p$.

We consider a mechanical oscillator with effective mass $m$ and natural
resonance frequency $\Omega_{0}$, weakly coupled to a thermal bath at
temperature $T$. This oscillator undergoes continuous position measurement,
and is concurrently manipulated by measurement-based feedback designed to
shift its resonance frequency from $\Omega_{0}$ to $\Omega$. The oscillator’s
position and momentum, $x$ and $p$, provide a description of its state. They
comply with the commutation relation $[x,p]=i\hbar$, fundamentally limiting
the precision with which they can be simultaneously known. When the oscillator
is cooled to its ground state, its position and momentum are localized to
their respective zero-point fluctuations,
$x_{\text{zp}}=\sqrt{{\hbar}/{2m\Omega}}$ and $p_{\text{zp}}=\sqrt{\hbar
m\Omega/2}$. For a squeezed state, one observable has uncertainty below the
zero-point fluctuations, with the uncertainty in the other correspondingly
increased.

To determine the uncertainty in the observables, we employ state estimation
following the approach in Ref. [8], but extending it to include feedback
control. The estimation relies on the measurement record, with the shared
measurement channel servicing both estimation and feedback control, as shown
in Fig. 1. The system features complex dynamics due to the interplay between
measurement conditioning, dissipation, bath noise, quantum backaction, and
measurement-based feedback force.

While our analysis is applicable more broadly, we focus on implementation
using an optomechanical system equipped with an optical cavity, boosting the
optomechanical interaction strength [16, 13], as shown in Fig. 1. We adhere to
the unresolved sideband regime, for which the optical cavity decay rate
$\kappa$ is much larger than both $\Omega_{0}$ and $\Omega$. In this limit,
both the mechanical oscillator and the input field can be adiabatically
tracked [17, 18]. The laser frequency is chosen to be on-resonance with the
cavity, causing the mechanical position to be imprinted on the output optical
phase quadrature [13, 16]. Measuring this quadrature provides an inference of
the mechanical position, $x_{\mathrm{meas}}(t)=x(t)+\delta
x_{\mathrm{imp}}(t)$, where $\delta x_{\mathrm{imp}}(t)$ is the imprecision
noise of the measurement. The optomechanical interaction also introduces
measurement backaction, imprinting the intracavity optical amplitude
quadrature onto the mechanical resonator. A feedback force proportional to
$x_{\mathrm{meas}}$ is applied to the mechanical oscillator. Experimentally,
this could be applied through modulations of the incident optical intensity
[19, 20, 21, 22], electrical actuation [23, 24, 25, 26, 27] or other methods.

We characterize the mechanical motion in the frequency domain by taking the
Fourier transform with convention
$F(\omega)=\int^{+\infty}_{-\infty}f(t)e^{i\omega t}dt$. The position response
is determined by the combined influence of the thermal force
$F_{\mathrm{th}}(\omega)$, backaction force $F_{\mathrm{ba}}(\omega)$, and
feedback force $F_{\mathrm{fb}}(\omega)$. This is given by [27]:

$x(\omega)=\left(F_{\mathrm{th}}(\omega)+F_{\mathrm{ba}}(\omega)+F_{\mathrm{fb}}(\omega)\right)\chi_{0}(\omega),$
(1)

where $\chi_{0}(\omega)$ is the intrinsic susceptibility defined by
$\left.\chi_{0}^{-1}(\omega)={m}{(\Omega_{0}^{2}-\omega^{2}-i\Gamma\omega)}\right.$
with $\Gamma$ being the mechanical energy damping rate. The feedback force,
which depends on the measurement, can be expressed as
$\left.F_{\mathrm{fb}}(\omega)=\chi_{\mathrm{fb}}^{-1}(\omega)x_{\mathrm{meas}}(\omega)\right.$.
In general, the feedback susceptibility $\chi_{\mathrm{fb}}(\omega)$ is a
complex function of frequency, taking into account any filtering of the
measured photocurrent, feedback strength and delays in the feedback loop
111$\chi_{\mathrm{fb}}^{-1}(\omega)=-Km(1-i\omega/\Gamma_{\mathrm{fb}})e^{i\omega\tau}\approx-
Km(1+\omega^{2}\tau/\Gamma_{\mathrm{fb}}+i\omega(\tau-1/\Gamma_{\mathrm{fb}}))$..
Here, for simplicity, we choose it to be a real constant,
$\left.\chi_{\mathrm{fb}}^{-1}(\omega)=-Km\right.$ with $K$ representing the
gain factor of the feedback circuit. The approximation holds for
$\omega\tau\ll 1$ [27], where $\tau$ is the overall delay in the feedback
loop. Choosing $\chi_{\mathrm{fb}}(\omega)$ to be real causes the feedback to
induce a mechanical resonance frequency shift. Solving for $x(\omega)$ under
this assumption, Eq. (1) can be re-written as
$(F_{\mathrm{th}}(\omega)+F_{\mathrm{ba}}(\omega)+\delta
F_{\mathrm{fb}}(\omega))\chi(\omega)=x(\omega)$, where $\delta
F_{\mathrm{fb}}(\omega)=\chi_{\mathrm{fb}}^{-1}(\omega)\delta
x_{\mathrm{imp}}(\omega)$ is a force noise term introduced by the measurement
imprecision. The modified mechanical susceptibility
$\left.\chi^{-1}(\omega)={m}{(\Omega^{2}-\omega^{2}-i\Gamma\omega)}\right.$,
with shifted mechanical frequency $\Omega=\sqrt{\Omega_{0}^{2}+K}$, and
unchanged mechanical dissipation rate.

We now estimate the mechanical state in the presence of both linear position
measurement and feedback control. The measurement signal satisfies the
commutation relation $[x_{\mathrm{meas}}(t),x_{\mathrm{meas}}(t^{\prime})]=0$
[29, 30]. This enables it to be treated classically and classical filtering to
be used to estimate the mechanical position and momentum. The estimates are
given by $\left.x_{\text{est}}(t)=H_{x}(t)\circledast
x_{\mathrm{mea}}(t)\right.$ and $\left.p_{\text{est}}(t)=H_{p}(t)\circledast
x_{\mathrm{mea}}(t)\right.$, where $H_{x}(t)$ and $H_{p}(t)$ are filter
functions. The optimal estimates, minimizing the mean-squared error, are
provided by the causal Wiener filters (see derivation in Ref.
100100100Supplemental material) [32]

${H}_{o}(\omega)=\left(A_{o}-iB_{o}\omega\right)\chi(\omega)^{\prime},\quad
o\in\\{x,p\\}$ (2)

where $A_{o}$ and $B_{o}$ are frequency-independent coefficients given in Ref.
100100footnotemark: 100, and $\left.\chi(\omega)^{\prime}={1}/{(\Omega^{\prime
2}-\omega^{2}-i\Gamma^{\prime}\omega)}\right.$ corresponds mathematically to a
modified mechanical susceptibility that peaks at
$\left.\Omega^{\prime}=\left(\Omega_{0}^{4}+\Omega_{\mathrm{meas}}^{4}\right)^{1/4}\right.$
and has linewidth
$\left.\Gamma^{\prime}=\sqrt{\Gamma^{2}-2\Omega_{0}^{2}+2\Omega^{\prime
2}}\right.$. The frequency $\Omega_{\mathrm{meas}}$ is an important parameter
that emerges from the model and is independent of feedback control. We term
this the characteristic measurement frequency. It is given by
$\left.\Omega_{\mathrm{meas}}=2(\eta
Cn_{\mathrm{tot}})^{1/4}\sqrt{\Gamma\Omega_{0}}\right.$, where $\eta$ is the
detection efficiency, $\left.n_{\mathrm{tot}}=n_{\mathrm{th}}+C+1/2\right.$ is
the intrinsic total occupancy (including the 1/2 quanta of vacuum energy for
succinctness), $n_{\mathrm{th}}\approx k_{B}T/\hbar\Omega$ is the average
thermal phonon occupation [33], and $C=4g^{2}/\kappa\Gamma$ is the
optomechanical cooperativity with $g$ the coherent amplitude boosted
optomechanical coupling rate of the un-frequency-shifted oscillator [13, 16].
For an oscillator without feedback control,
$\left.\Omega_{\mathrm{meas}}>\Omega_{0}\right.$ defines the regime in which
the rotating wave approximation (RWA) is invalid [8].

As the measurement is linear, the oscillator’s conditional state can be fully
described by its covariance matrix,
$\mathbb{V}=\left(\begin{array}[]{cc}V_{\delta X\delta X}&V_{\delta X\delta
P}\\\ V_{\delta X\delta P}&V_{\delta P\delta P}\end{array}\right),$ where here
and throughout we normalize the variances and covariances by their respective
zero-point motion (i.e., $\left.V_{\delta X\delta X}=V_{\delta x\delta
x}/x_{\mathrm{zp}}^{2}\right.$, $\left.V_{\delta P\delta P}=V_{\delta p\delta
p}/p_{\mathrm{zp}}^{2}\right.$, $\left.V_{\delta X\delta P}=V_{\delta x\delta
p}/x_{\mathrm{zp}}p_{\mathrm{zp}}\right.$, and $\left.{\delta x,\delta
p}={x-x_{\mathrm{est}},p-p_{\mathrm{est}}}\right.$). Employing the optimal
Wiener filters, we derive analytical expressions for the covariance matrix.
These are characterized by five dimensionless variables: $\eta$, $C$,
$n_{\mathrm{th}}$, $Q$, and the relative frequency ratio
$R=\Omega/\Omega_{0}$. In the general case, the expressions are non-intuitive.
For completeness, we include them in Ref. 100100footnotemark: 100.

Fig. 2 plots the diagonal terms of the covariance matrix as a function of $C$
for different choices of feedback gain. In the absence of feedback control,
both mechanical position and momentum variances decrease as $C$ increases in
the thermal noise dominated regime ($C\ll n_{\mathrm{th}}$). The variances
stabilize, becoming equally localized close to the zero-point fluctuations in
the backaction-dominated regime ($C\gg n_{\mathrm{th}}$). Further increasing
$C$ into the squeezing regime ($C>Q/2$ 222Inside the backaction dominated
regime, the squeezing criterion $C>\left(n_{\mathrm{tot}}Q^{2}/4\right)^{1/3}$
can be simplified to $\left.C>Q/2\right.$) leads to a reduction in position
variance below the zero-point fluctuations, while the momentum variance
increases, consistent with predictions in previous work [8].

Feedback control can influence the conditional state without altering the
intrinsic backaction and thermal noise. As shown in Fig. 2, in the softening
scenario ($R<1$) we find that the position localization is enhanced and the
momentum localization is degraded. As $R$ decreases, the conditional position
variance can be reduced beneath the zero-point fluctuations even in scenarios
where the no-feedback quantum squeezing criterion defined in Ref. [8] is not
satisfied. In complement to this, the hardening scenario ($R>1$) weakens
position localization, making quantum squeezing more challenging to achieve.
Surprisingly, given the use of position measurement, this hardening allows
momentum squeezing. This is noteworthy given that momentum rather than
position squeezing is required for quantum enhanced force sensing [35]. The
conditional momentum variance initially decreases with increasing $C$ and
eventually increases for sufficient high $C$. Thus, unlike position squeezing,
which monotonically improves as $C$ increases, we find that there is an
optimal $C$ for maximal momentum squeezing.

Figure 2:  Characterization of conditional states under varying feedback with
$n_{\mathrm{th}}=2\times 10^{3}$ and $\eta=1$ as a function of $C$. (a)
Position and (b) momentum variances for several $R$, are shown with cold and
warm colors, respectively. Grey area: squeezing regime. Red dashed lines:
backaction dominated ($C>n_{\mathrm{th}}$) and quantum squeezing boundaries
($C>Q/2$) without feedback.

To provide a comprehensive understanding of how the feedback control affects
the squeezed state in different measurement regimes, we characterize the
conditional position and momentum variances as a function of $n_{\mathrm{th}}$
and $C$ in Fig. 3.

Fig. 3 (a) illustrates the position variance within the softening scenario for
$R=0.1$. Importantly, the parameter regime in which quantum squeezing occurs
significantly expands when feedback is applied (blue c.f. grey shaded
regions), enabling squeezing in regimes previously deemed inaccessible. This
eases the requirements for achieving quantum squeezing, opening the
possibility for quantum state preparation with a significantly wider class of
optomechanical systems.

Fig. 3 (b) presents position (blue) and momentum (orange) variances within the
hardening regime for $R=10$. Notably, the parameter regime supporting position
squeezing contracts compared to the case without feedback. However, this
contraction is accompanied by the emergence of a new regime characterized by
momentum squeezing as shown in the orange region in Fig. 3 (b).

To further understand the impact of feedback on the conditional state and its
subsequent influence on the criteria for quantum squeezing, we simplify the
expressions governing the conditional state by taking the limit of a high $Q$
oscillator ($Q\gg 1$) and finite measurement strength, i.e., $C\not\to\ 0$.
The covariance matrix of the conditional state can then be greatly simplified
to

$\mathbb{V}=\frac{2n_{\text{tot }}\Gamma}{\Omega_{\text{meas
}}^{4}}\left(\begin{array}[]{cc}\Gamma^{\prime}R\Omega_{0}^{2}&{2\Omega_{0}}\left(\Omega^{\prime
2}-\Omega_{0}^{2}\right)\\\ {2\Omega_{0}}\left(\Omega^{\prime
2}-\Omega_{0}^{2}\right)&\Gamma^{\prime}\Omega^{\prime
2}/R\end{array}\right).$ (3)

The purity of the system $\mathcal{P}=\sqrt{{\eta C}/{n_{\text{tot}}}}$ is
independent of feedback control and is always equal to or below unity,
ensuring the state’s physical validity. The diagonal terms in Eq. (3)
characterize the position and momentum squeezing. The non-zero off-diagonal
covariance terms indicate that the position and momentum are correlated so
that the optimal squeezing occurs at an in-between quadrature (detailed in
Ref. 100100footnotemark: 100). However, for typical measurement regimes, the
variance of this quadrature is close to the minimum of $V_{\delta x\delta x}$
and $V_{\delta p\delta p}$. We thus confine our analysis to position and
momentum.

Figure 3:  Characterization of the conditional variances as a function of
$n_{\mathrm{th}}$ and $C$ for $R=0.1$ (a) and $R=10$ (b), with $\eta=1$, at
room temperature. Black dashed lines: RWA breakdown criterion. Dotted lines in
corresponding colors: respective squeezing thresholds. Grey shaded regions:
position squeezing regimes when $\left.R=1\right.$. Blue shading: position
squeezing. Orange shading: momentum squeezing.

We establish intuitive criteria for position and momentum squeezing by
simplifying Eq. (3) and determining when its diagonal elements are less than
unity. To do this, we consider two distinct scenarios: the non-RWA regime
$\Omega_{\text{meas}}\gg\Omega_{0}$, for which
$\Omega^{\prime}\rightarrow\Omega_{\text{meas}}$ and
$\Gamma^{\prime}\rightarrow\sqrt{2}\Omega_{\text{meas}}$; and the RWA regime
$\Omega_{\text{meas}}\ll\Omega_{0}$, for which
$\Omega^{\prime}\rightarrow\Omega_{0}$ and
$\Gamma^{\prime}\rightarrow\Omega_{\text{meas}}^{2}/\Omega$. To facilitate
comparison with optomechanical systems that operate without feedback control,
we formulate these criteria using unshifted dimensionless parameters $\eta$,
$C$, $n_{\mathrm{th}}$ and $Q$, consistent with extant literature [8].

In the non-RWA regime, we find that position squeezing is achieved when
$C>\left(n_{\mathrm{tot}}Q^{2}R^{4}/4\right)^{1/3}/\eta$. This criterion can
be further simplified if either thermal noise ($n_{\mathrm{th}}\gg C$) or
backaction noise ($C\gg n_{\mathrm{th}}$) dominates the noise driving the
mechanical oscillator. In the former case, we find that squeezing occurs when
$\left.C>\left(n_{\mathrm{th}}Q^{2}R^{4}/4\right)^{1/3}/\eta\right.$, while in
the latter it occurs when $\left.C>QR^{2}/2\eta^{3/2}\right.$. Hence, feedback
reduces the required $C$ for squeezing by factors of $R^{4/3}$ and $R^{2}$,
respectively. In an ideal scenario with perfect detection efficiency
($\left.\eta=1\right.$) and in the backaction dominated regime, the position
squeezing criterion reduces to
$\left.\Omega_{\mathrm{meas}}>\sqrt{2}R\Omega_{0}\right.$, so that position
squeezing is achievable when the measurement bandwidth surpasses the feedback-
shifted mechanical frequency by a factor of $\sqrt{2}$.

Within the RWA regime, our findings reveal that squeezing can be achieved when
$C>n_{\mathrm{tot}}R^{2}/\eta$, independent of whether the system is in the
backaction or thermal noise dominated regimes.

Analogous to position squeezing, we establish the momentum squeezing criteria
as $C>64n_{\mathrm{tot}}^{3}/\eta R^{4}Q^{2}$ and $C>n_{\mathrm{tot}}/\eta
R^{2}$ outside and inside the RWA regimes, respectively. Notably, feedback not
only makes momentum squeezing possible, but can also be used to tune the
requirements for squeezing through adjustable feedback strength.

The position and momentum squeezing criteria, deduced in the previous
paragraphs are depicted by blue and orange dotted lines in Fig. 3 panels (a)
and (b), agreeing well with the full model.

Figure 4: Conditional position variance as a function of $n_{\mathrm{cav}}$
for different $R$ at room temperature, with (colored lines) and without (black
line) second-order mode noise. The split-colored lines represent the lower and
upper bounds of the conditional variance when taking the simplified model.
Fundamental mode optomechanical parameters: $G/2\pi\sim$14 MHz/nm,
$\Omega_{0}/2\pi\sim$0.8 MHz, $Q\sim~{}10^{6}$ and $m\sim$7 ng. Second mode:
$G_{2}=G$, $\Omega_{2}=2\Omega_{0}$, $Q_{2}=Q$, and $m_{2}=m$. $\eta$=0.63
[36].

To evaluate the experimental feasibility of our proposed method, we consider
the membrane-in-the-middle optomechanical system of Ref. [36]. This square
membrane has optomechanical coupling strength $G/2\pi\sim$14 MHz/nm, detection
efficiency $\eta$=0.63, and a fundamental mode with $\Omega_{0}/2\pi\sim$0.8
MHz, $Q\sim~{}10^{6}$ and $m\sim$7 ng [36]. With these parameters, at room
temperature position squeezing can be achieved with $3\times 10^{9}$
intracavity photons. With feedback to achieve $R=1/10$, we find this to reduce
to $9\times 10^{7}$ intracavity photons. By contrast, increasing the resonance
frequency by a factor of ten does not enable momentum squeezing. Momentum
squeezing is possible though, for higher $R$. For $R=20$, it can be achieved
with $3\times 10^{8}$ intracavity photons. Ref. [36] used $4\times 10^{8}$
intracavity photons, indicating that feedback could realistically enable both
position and momentum squeezing.

Background mechanical resonances have been identified to significantly degrade
measurement-based state preparation in a recent experiment [15]. It is
interesting to examine if feedback mechanisms can counteract this effect and
facilitate mechanical squeezing. To explore this, we examine the influence of
the second higher-order mode of the membrane considered above. We use the
parameters $\Omega_{2}=2\Omega_{0}$, $Q_{2}=Q$, $m_{2}=m$, and a conservative
estimate for optomechanical coupling strength of $G_{2}=G$ [36, 37]. The noise
of the second resonance peak has a flat spectrum at frequencies well below
$\Omega_{2}$ and is suppressed by the Wiener filter at frequencies well above
$\Omega_{0}$ [15]. As such, it can be well-approximated as an additional white
noise source. This simplifies the calculation of the conditional variance,
allowing analytical solutions. The amplitude of the additional noise is
bounded from below by the noise from the second mode at zero frequency, and
from above by the maximum noise from the second mode across the full frequency
band that the Wiener filter samples. This results in lower and upper bounds,
respectively, on the conditional variance, as further justified in Ref.
100100footnotemark: 100.

Contrary to the case of a single mechanical mode, for which the position
variance decreases monotonically with measurement strength [8], within the
validity of the white noise approximation our analysis shows that the presence
of a second mode introduces an optimal measurement strength, above which the
variance degrades. This is shown in Fig. 4. The optimum measurement is
independent of feedback, and for the chosen parameters, this optimum occurs
around $10^{9}$ intracavity photons. As shown in Fig. 4, even at the optimum
intracavity photon number, squeezing is not possible in the absence of
feedback. However, it is possible with feedback for $R\lesssim 1/25$.

In conclusion, we have shown that feedback control can substantially relax the
requirements for measurement-based mechanical position squeezing, and that it
can enable momentum squeezing. Furthermore, we find that feedback control can
suppress degradation in squeezing due to the presence of background mechanical
resonances. This paves the way for room-temperature quantum state preparation
and sensing.

###### Acknowledgements.

We would like to thank J. X. Chen, S. A. Fedorov, J. S. Bennett, A.
Schliesser, K. Mølmer, E. Zeuthen, E. S. Polzik for useful discussion. This
work was supported by the Australian Research Council Centres of Excellence
for Engineered Quantum Systems (EQUS, CE170100009) and Quantum Biotechnology
(CE230100021), Air Force Office of Scientific Research under award number
FA9550-20-1-0391 and VILLUM FONDEN under a Villum Investigator Grant No.
25880. C.M. acknowledges funding from the European Union’s Horizon Europe
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No. 101110196.

Supplemental Material: Enhancement of mechanical squeezing via feedback
control

##  I Model of measurement and feedback

Estimating the state utilizing optimal Wiener filters necessitates a
comprehensive system model. This model is effectively represented through the
spectral densities of displacement and force noise. The displacement spectral
density is formulated as [27]

$S_{xx}(\omega)=(S_{FF}^{\text{tot }}(\omega)+S_{FF}^{\mathrm{\delta
fb}}(\omega))|\chi(\omega)|^{2},$ (S1)

where $S_{FF}^{\mathrm{tot}}(\omega)$ encapsulates the total spectral density
of both thermal $S_{FF}^{\mathrm{th}}(\omega)$ and backaction
$S_{FF}^{\mathrm{ba}}(\omega)$ forces. This is expressed as

$S_{FF}^{\mathrm{tot}}(\omega)=2n_{\mathrm{tot}}S_{xx}^{\mathrm{zp}}(\Omega_{0})/|\chi_{0}(\Omega_{0})|^{2},$
(S2)

where
$S_{xx}^{\mathrm{zp}}(\Omega_{0})=2x_{\mathrm{zp}}(\Omega_{0})^{2}/\Gamma$
defines the intrinsic zero-point noise of the oscillator [38, 27]. The
spectral density characterizing the imprecision feedback force given by
$S_{FF}^{\mathrm{\delta
fb}}(\omega)=S_{xx}^{\mathrm{imp}}(\omega)|\chi_{\mathrm{fb}}(\omega)|^{-2}$,
where
$\left.S_{xx}^{\mathrm{imp}}(\omega)={S_{xx}^{\mathrm{zp}}(\Omega_{0})}/{8\eta
C}\right.$ indicates the measurement-imprecision noise at the optical vacuum
noise level.

Subsequently, the noise equivalent Power Spectral Density (PSD) of the
measurement laser phase quadrature is given by [38, 16, 27]

$S_{xx}^{\mathrm{meas}}(\omega)=S_{xx}(\omega)+S_{xx}^{\mathrm{imp}}(\omega)+2S_{x\delta
x_{\mathrm{imp}}}(\omega),$ (S3)

where $S_{x\delta x_{\mathrm{imp}}}(\omega)$ represents the cross-correlation
between the displacement and imprecision noise introduced by the feedback. The
cross-spectral density is defined by [38, 16]

$S_{AB}(\omega)=\int^{\infty}_{-\infty}e^{i\omega t}\langle
A(t)B(0)\rangle\text{d}t=\int^{\infty}_{-\infty}\frac{d\omega^{\prime}}{2\pi}\langle{A^{{\dagger}}(-\omega)}{B(\omega\prime)}\rangle.$

##  II Measurement bandwidth

In the high $Q$ limit where $Q=\Omega_{0}/\Gamma\gg 1$, $\Omega_{\text{meas}}$
can be directly related to the measurement bandwidth $B_{\text{meas}}$ for
which the intrinsic thermal noise is resolved above the optical shot noise.
Inside the RWA regime,
$\Omega_{\mathrm{meas}}=\sqrt{B_{\text{meas}}\Omega_{0}}$, while outside this
regime $\left.\Omega_{\mathrm{meas}}=B_{\text{meas}}\right.$.

##  III Filter functions

The causal Wiener filter is [32]

$H_{o}(\omega)=\frac{1}{M_{x}^{\mathrm{meas}}(\omega)}\left[\frac{S_{ox^{\mathrm{meas}}}(\omega)}{{M_{x}^{\mathrm{meas}}(\omega)^{*}}}\right]_{+},\quad
o\in\\{x,p\\}$ (S4)

where $S_{ox^{\mathrm{meas}}}(\omega)=S_{ox}(\omega)+S_{o\delta
x_{\mathrm{imp}}}(\omega)$, and $M_{x}^{\mathrm{meas}}(\omega)$ is the causal
spectral factor that satisfies
$\left.S_{xx}^{\mathrm{meas}}(\omega)=M_{x}^{\mathrm{meas}}(\omega)M_{x}^{\mathrm{meas}(\omega)^{*}}\right.$
and only has poles and zeros in the lower half of the complex plane.
$[...]_{+}$ denotes the causal part of the function. In general, any function
may be separated into the sum of its causal and anti-causal parts [32], which
can be decomposed by factorizing the poles of the denominator into the upper
and the lower halves of the complex plane and finding the partial fraction
decomposition. Expanding Eq. (S4) using this procedure, we find the
coefficients of the filter functions in Eq. (2) are

$\displaystyle A_{x}=\frac{\Omega_{\text{meas
}}^{4}\left(\Gamma^{2}+\Gamma\Gamma^{\prime}+\Omega^{\prime
2}-\Omega^{2}\right)+\left(\Omega_{0}^{2}-\Omega^{2}\right)\left(\Omega_{0}^{2}\left(\Gamma\left(\Gamma+\Gamma^{\prime}\right)+\Omega^{\prime
2}\right)-\left(\Omega_{0}^{2}+\Omega^{\prime
2}\right)\Omega^{2}+\Omega^{4}\right)}{\left(\Gamma\left(\Gamma+\Gamma^{\prime}\right)\Omega^{\prime
2}+\Omega^{\prime
4}+\left(\Gamma^{\prime}\left(\Gamma+\Gamma^{\prime}\right)-2\Omega^{\prime
2}\right)\Omega^{2}+\Omega^{4}\right)},$ (S5) $\displaystyle
B_{x}=\frac{\Omega_{\text{meas
}}^{4}\left(\Gamma+\Gamma^{\prime}\right)+\left(\Omega_{0}^{2}-\Omega^{2}\right)\left(\Gamma\left(\Omega_{0}-\Omega^{\prime}\right)\left(\Omega_{0}+\Omega^{\prime}\right)+\Gamma^{\prime}\left(\Omega_{0}^{2}-\Omega^{2}\right)\right)}{\left(\Gamma\left(\Gamma+\Gamma^{\prime}\right)\Omega^{\prime
2}+\Omega^{\prime
4}+\left(\Gamma^{\prime}\left(\Gamma+\Gamma^{\prime}\right)-2\Omega^{\prime
2}\right)\Omega^{2}+\Omega^{4}\right)},$ $\displaystyle
A_{p}=m\frac{\Omega^{2}\left(-\Omega_{\text{meas
}}^{4}\left(\Gamma+\Gamma^{\prime}\right)-\left(\Omega_{0}^{2}-\Omega^{2}\right)\left(\Gamma\left(\Omega_{0}^{2}-\Omega^{\prime
2}\right)+\Gamma^{\prime}\left(\Omega_{0}^{2}-\Omega^{2}\right)\right)\right)}{\left(\Gamma\left(\Gamma+\Gamma^{\prime}\right)\Omega^{\prime
2}+\Omega^{\prime
4}+\left(\Gamma^{\prime}\left(\Gamma+\Gamma^{\prime}\right)-2\Omega^{\prime
2}\right)\Omega^{2}+\Omega^{4}\right)},$ $\displaystyle
B_{p}=m\frac{\Omega_{\text{meas }}^{4}\left(\Omega^{\prime
2}-\Omega^{2}\right)+\left(\Omega_{0}^{2}-\Omega^{2}\right)\left(\Gamma^{2}\Omega^{\prime
2}+\Gamma\Gamma^{\prime}\Omega^{2}+\left(\Omega^{\prime
2}-\Omega^{2}\right)\left(\Omega_{0}^{2}-\Omega^{2}\right)\right)}{\left(\Gamma\left(\Gamma+\Gamma^{\prime}\right)\Omega^{\prime
2}+\Omega^{\prime
4}+\left(\Gamma^{\prime}\left(\Gamma+\Gamma^{\prime}\right)-2\Omega^{\prime
2}\right)\Omega^{2}+\Omega^{4}\right)}.$

##  IV Covariance matrix

Employing the optimal filters in Eq. (2), we derive the variances and
covariance using [16]

$\displaystyle V_{\delta x\delta
x}=\int^{\infty}_{-\infty}\frac{d\omega}{2\pi}S_{\delta x\delta x},$ (S6)
$\displaystyle V_{\delta p\delta
p}=\int^{\infty}_{-\infty}\frac{d\omega}{2\pi}S_{\delta p\delta p},$
$\displaystyle V_{\delta x\delta
p}=\int^{\infty}_{-\infty}\frac{d\omega}{2\pi}\mathrm{Re}\\{S_{\delta x\delta
p}\\}.$

The expressions for the normalized, fully analytical variances and covariance
of conditional position and momentum are

$\displaystyle V_{\delta X\delta
X}=\frac{1}{8CR\Gamma^{2}}\left(16Cn_{\text{tot
}}\Gamma^{2}+\left(-1+R^{2}\right)^{2}\Omega_{0}^{2}+\left(\Omega_{0}^{2}\left(-2R^{2}\left(-1+R^{2}\right)\Omega_{0}^{2}\left(16Cn_{\text{tot
}}\Gamma^{2}-\left(-1+R^{2}\right)\Omega_{0}^{2}\right)\left(R^{4}\Omega_{0}^{4}+\right.\right.\right.\right.$
(S7)
$\displaystyle\left.2\Gamma\left(\Gamma+\Gamma^{\prime}\right)\Omega^{\prime
2}+\Omega^{\prime 4}+R^{2}\Omega_{0}^{2}\left(-\Gamma^{2}+\Gamma^{\prime
2}-2\Omega^{\prime
2}\right)\right)-R^{2}\left(-1+R^{2}\right)^{2}\Omega_{0}^{2}(R^{6}\Omega_{0}^{6}+\Gamma^{2}\Omega^{\prime
4}+R^{4}\Omega_{0}^{4}\left(\Gamma^{\prime 2}-2\Omega^{\prime 2}\right)+$
$\displaystyle R^{2}\Omega_{0}^{2}\Omega^{\prime 2}\left(2\Gamma\Gamma
1+\Omega^{\prime 2}\right))-\left(-16Cn_{\text{tot
}}\Gamma^{2}+\left(-1+R^{2}\right)\Omega_{0}^{2}\right)^{2}\left.\left.\left(R^{4}\Omega_{0}^{4}+R^{2}\Omega_{0}^{2}\left(-\Gamma^{2}+\Gamma^{\prime
2}-2\Omega^{\prime
2}\right)+\left(\Gamma\left(\Gamma+\Gamma^{\prime}\right)+\Omega^{\prime
2}\right)^{2}\right)\right)\right)/$
$\displaystyle\left.\left(R^{4}\Omega_{0}^{4}+\Gamma\left(\Gamma+\Gamma^{\prime}\right)\Omega^{\prime
2}+\Omega^{\prime 4}+R^{2}\Omega_{0}^{2}\left(\Gamma
1^{\prime}\left(\Gamma+\Gamma^{\prime}\right)-2\Omega^{\prime
2}\right)\right)^{2}\right)$ $\displaystyle V_{\delta P\delta
P}=\left(16Cn_{\text{tot}}\Gamma^{2}\left(-2R^{2}\left(\Gamma+\Gamma^{\prime}\right)^{2}\Omega_{0}^{6}-R^{8}\Omega_{0}^{8}+\right.\right.R^{4}\Omega_{0}^{4}\left(\Gamma^{\prime
2}\left(\Gamma+\Gamma^{\prime}\right)^{2}+2\left(\Gamma^{2}+4\Gamma\Gamma^{\prime}+2\Gamma^{\prime
2}\right)\Omega_{0}^{2}-2\Omega_{0}^{4}\right)+R^{6}(-2\Gamma\Gamma^{\prime}\Omega_{0}^{6}+$
$\displaystyle
4\Omega_{0}^{8})+2R^{2}\Omega_{0}^{2}\left(\Gamma\Gamma^{\prime}\left(\Gamma+\Gamma^{\prime}\right)^{2}-\left(\left(-2+R^{2}\right)\Gamma^{2}+R^{2}\Gamma\Gamma^{\prime}+2R^{2}\Gamma^{\prime
2}\right)\Omega_{0}^{2}+2\left(1-2R^{2}\right)\Omega_{0}^{4}\right)\Omega^{\prime
2}+(\Gamma^{2}\left(\Gamma+\Gamma^{\prime}\right)^{2}-2(\left(1+R^{2}\right)\Gamma^{2}$
$\displaystyle+R^{2}\Gamma\Gamma^{\prime}-R^{2}\Gamma^{\prime
2})\Omega_{0}^{2}+2\left(-1+2\left(R^{2}+R^{4}\right)\right)\Omega_{0}^{4})\Omega^{\prime
4}+\left.2\left(\Gamma\left(\Gamma+\Gamma^{\prime}\right)-2R^{2}\Omega_{0}^{2}\right)\Omega^{\prime
6}+\Omega^{\prime 8}\right)-$ $\displaystyle 256C^{2}n_{\text{tot
}}^{2}\Gamma^{4}\Omega_{0}^{2}\left(R^{4}\Omega_{0}^{4}+\Omega^{\prime
4}+R^{2}\Omega_{0}^{2}\left(\left(\Gamma+\Gamma^{\prime}\right)^{2}-2\Omega^{\prime
2}\right)\right)+\left(-1+R^{2}\right)^{2}\Omega_{0}^{2}\left(\Omega^{\prime
4}\left(\Gamma\Gamma^{\prime}-\Omega_{0}^{2}+\Omega^{\prime
2}\right)\left(\Gamma\left(2\Gamma+\Gamma^{\prime}\right)+\Omega_{0}^{2}+\Omega^{\prime
2}\right)+\right.$ $\displaystyle
R^{2}\Omega_{0}^{2}\left(-\left(\Gamma+\Gamma^{\prime}\right)^{2}\Omega_{0}^{4}+2\left(\Gamma\Gamma^{\prime
3}+\Omega_{0}^{4}+2\Gamma^{2}\left(\Gamma^{\prime
2}+\Omega_{0}^{2}\right)\right)\Omega^{\prime
2}+\right.\left.\left(-3\Gamma^{2}-2\Gamma\Gamma^{\prime}+2\left(\Gamma^{\prime
2}+\Omega_{0}^{2}\right)\right)\Omega^{\prime 4}-4\Omega^{\prime 6}\right)$
$\displaystyle+R^{6}\Omega_{0}^{6}\left(\Gamma^{\prime
2}+2\left(\Omega_{0}-\Omega^{\prime}\right)\left(\Omega_{0}+\Omega^{\prime}\right)\right)+\left.\left.R^{4}\Omega_{0}^{4}\left(\Gamma^{\prime
4}+2\Gamma^{\prime 2}\Omega_{0}^{2}-\Omega_{0}^{4}-4\left(\Gamma^{\prime
2}+\Omega_{0}^{2}\right)\Omega^{\prime 2}+5\Omega^{\prime
4}+2\Gamma\Gamma^{\prime}\left(\Gamma^{\prime
2}+2\Omega_{0}^{2}-\Omega^{\prime 2}\right)\right)\right)\right)/$
$\displaystyle\left(8CR\Gamma^{2}\left(R^{4}\Omega_{0}^{4}+\Gamma\left(\Gamma+\Gamma^{\prime}\right)\Omega^{\prime
2}+\Omega^{\prime
4}+R^{2}\Omega_{0}^{2}\left(\Gamma^{\prime}\left(\Gamma+\Gamma^{\prime}\right)-2\Omega^{\prime
2}\right)\right)^{2}\right)$ $\displaystyle V_{\delta X\delta
P}=\frac{\Omega_{0}^{3}\left(16Cn_{\text{tot
}}\Gamma^{2}\left(\Gamma+\Gamma^{\prime}\right)+(-1+R)(1+R)\left(-\left(\left(\Gamma+\Gamma^{\prime}-R^{2}\Gamma^{\prime}\right)\Omega_{0}^{2}\right)+\Gamma\Omega^{\prime
2}\right)\right)^{2}}{8C\Gamma\left(R^{4}\Omega_{0}^{4}+\Gamma\left(\Gamma+\Gamma^{\prime}\right)\Omega^{\prime
2}+\Omega^{\prime
4}+R^{2}\Omega_{0}^{2}\left(\Gamma^{\prime}\left(\Gamma+\Gamma^{\prime}\right)-2\Omega^{\prime
2}\right)\right)^{2}}.$

In the main text of our analysis, Eq. (S7) has been simplified to Eq. (3) by
considering a high $Q$ oscillator. Further simplification of Eq. (3) is
achievable, both within and outside the RWA regime.

Far outside the RWA regime ($\Omega_{\mathrm{meas}}\gg\Omega_{0}$), the
covariance matrix derived in Eq. (3) can then be simplified to

$\mathbb{V}_{\Omega_{\mathrm{meas}}\gg\Omega_{0}}=\frac{1}{\mathcal{P}}\left(\begin{array}[]{cc}\sqrt{2}\Omega/\Omega_{\text{meas
}}&1\\\ 1&\sqrt{2}\Omega_{\text{meas }}/\Omega\end{array}\right).$ (S8)

Conversely, deeply inside the RWA regime
($\Omega_{\mathrm{meas}}\ll\Omega_{0}$), yet outside the weak measurement
regime 333The protocol does not accommodate the weak measurement regime, as
the thermal motion cannot be resolved from the measurement, thereby rendering
feedback application unfeasible. Furthermore, it is worth noting that, through
feedback control, shifting $\Omega$ down nearly to DC only requires
$\Omega_{\mathrm{meas}}\gg\sqrt{\Omega_{0}\Gamma}$ while ground state cooling
requires $\Omega_{\mathrm{meas}}\gg
2n_{\mathrm{th}}n_{\mathrm{tot}}\sqrt{\Omega_{0}\Gamma}$, i.e., $\eta C\gg
n_{\mathrm{th}}$. Therefore, feedback softening does not have as significant a
technical hurdle as feedback cooling.
($\left.\Omega_{\mathrm{meas}}\gg\sqrt{\Omega_{0}\Gamma}\right.$) [16], the
correlation is suppressed so that in this regime $\theta=0$ and the covariance
matrix simplifies to

$\mathbb{V}_{\Omega_{\mathrm{meas}}\ll\Omega}=\frac{1}{\mathcal{P}}\left(\begin{array}[]{cc}R&0\\\
0&R^{-1}\end{array}\right).$ (S9)

## Minimum variance

The minimum variance $V_{\delta X^{\theta}\delta X^{\theta}}$ of the optimal
quadrature $\left.\delta X^{\theta}=\delta X\mathrm{cos}\theta+\delta
P\mathrm{sin}\theta\right.$ occurs for an angle
$\left.\theta=-\mathrm{arctan}(\sqrt{2}/(\Omega_{\mathrm{meas}}/\Omega-\Omega/\Omega_{\mathrm{meas}}))/2\right.$.
This minimum variance can be expressed in terms of the covariance matrix as
[8]

$V_{\delta X^{\theta}\delta X^{\theta}}=\frac{V_{\delta X\delta X}+V_{\delta
P\delta P}-\sqrt{V_{\delta X\delta X}^{2}+V_{\delta P\delta P}^{2}+4V_{\delta
X\delta P}^{2}-2V_{\delta X\delta X}V_{\delta P\delta P}}}{2}.$ (S10)

By substituting Eq. (S8) into Eq. (S10), we derive the simplified form

$V_{\delta X^{\theta}\delta X^{\theta}}=\frac{V_{\delta X\delta X}+V_{\delta
P\delta P}-\sqrt{V_{\delta X\delta X}^{2}+V_{\delta P\delta P}^{2}}}{2}.$
(S11)

In the regime of a strongly squeezed state, where
$\left.\mathrm{min}\\{V_{\delta X\delta X},V_{\delta P\delta P}\\}\ll
1\right.$, a first-order Taylor expansion applied to Eq. (S11) yields:

$V_{\delta X^{\theta}\delta X^{\theta}}=\frac{\mathrm{min}\\{V_{\delta X\delta
X},V_{\delta P\delta P}\\}}{2}.$ (S12)

##  V Model of higher order mode

In state preparation, higher-order mode noise detrimentally affects the
conditional state [15]. We model this effect by treating the higher-order mode
as white noise. To understand this, we use the error spectrum of the non-
causal Wiener filtering error $\Delta x=x-x_{\mathrm{est}}^{\mathrm{non}}$,
which has the feature $\Delta x\sim\delta x$ and is given by [32, 16, 22]

$S_{\Delta x\Delta x}(\omega)=\frac{S_{xx}(\omega)}{1+\mathrm{SNR}(\omega)},$
(S13)

with the signal-to-noise ratio $\mathrm{SNR}(\omega)$ defined as
$S_{xx}(\omega)/S_{nn}^{\mathrm{tot}}(\omega)$, where
$S_{nn}^{\mathrm{tot}}=S_{xx}^{\mathrm{imp}}(\omega)+S_{xx,2}(\omega)$ and
$S_{xx,2}(\omega)$ represents the PSD of the second-order mode noise. For high
and low $\mathrm{SNR}(\omega)$ regimes, Eq. (S13) simplifies to

$S_{\Delta x\Delta
x}(\omega)\approx\begin{cases}&S_{nn}^{\mathrm{tot}}(\omega)\quad\text{if}\quad\mathrm{SNR}(\omega)\gg
1,\\\ &S_{xx}(\omega)\quad\text{if}\quad\mathrm{SNR}(\omega)\ll 1.\end{cases}$
(S14)

As a result, we can reexpress

$S_{\Delta x\Delta
x}(\omega)\approx\mathrm{min}\\{S_{xx}(\omega),S_{nn}^{\mathrm{tot}}(\omega)\\}.$
(S15)

The error spectrum is the minimum of either the signal from the fundamental
mode or the total noise at a given frequency. When the intracavity photon
number is low that the condition $S_{xx,2}(\Omega_{12})<S_{xx}^{\mathrm{imp}}$
holds, where $\Omega_{12}$ represents the frequency at which the noise levels
of the fundamental and second-order modes intersect, i.e.,
$\left.S_{xx}(\Omega_{12})=S_{xx,2}(\Omega_{12})\right.$. In this scenario,
considering the noise from the second-order mode as white noise does not
affect the error spectrum, as shown in Fig. S1 (a). However, as the
intracavity photon number increases and reaches the regime where
$S_{xx}(\Omega_{12})>S_{xx}^{\mathrm{imp}}$, treating the second-order mode
noise as white at the magnitudes of $S_{xx,2}(0)$ and $S_{xx,2}(\Omega_{12})$
establishes the lower and upper bounds of the error spectrum, as illustrated
in Fig. S1 (b) and (c), respectively.

Figure S1:  Noise equivalent power spectral density using the square membrane
parameters. Purple line: signal from the fundamental mode. Green line: noise
from the second-order mode. Darkblue line: noise from the second-order mode at
its zero frequency noise level. Lightblue line: noise from the second-order
mode at its $\Omega_{12}$ frequency noise level. Grey line: optical shot
noise. The optomechanical parameters are chosen as same as in the main text.
(a) $S_{xx}(\Omega_{12})<S_{xx}^{\mathrm{imp}}$ for $n_{\mathrm{cav}}=10^{7}$.
Pink shaded area: $\mathrm{min}\\{S_{xx},S_{nn}^{\mathrm{tot}}\\}$. Diagonal
line pattern: the complete overlap of
$\mathrm{min}\\{S_{x}x,S_{xx}^{\mathrm{imp}}+S_{xx}(\Omega_{12})\\}$ and
$\mathrm{min}\\{S_{xx},S_{xx}^{\mathrm{imp}}+S_{xx}(0)\\}$. (b) and (c)
$S_{xx}(\omega_{12})>S_{xx}^{\mathrm{imp}}$ for $n_{\mathrm{cav}}=10^{9}$.
Diagonal line pattern:
$\mathrm{min}\\{S_{xx},S_{xx}^{\mathrm{imp}}+S_{xx}(0)\\}$ for (b) and
$\mathrm{min}\\{S_{xx},S_{xx}^{\mathrm{imp}}+S_{xx}(\Omega_{12})\\}$ for (c).

## Radiation pressure force and mechanical spring force

In the main text, one of the feedback control mechanisms we consider is the
radiation pressure force [19, 20, 21, 22]. It is worth considering whether the
radiation pressure force can be sufficiently strong to appreciably soften the
oscillator. To evaluate this, given the measurement strength, we compare the
radiation pressure force against the mechanical spring force. The latter
arises from thermal and backaction fluctuations. These forces are quantified
as $\left.F_{\text{rad }}=n_{\text{cav }}\hbar G\right.$ for the radiation
pressure and $\left.F_{\text{mech
}}=m\Omega_{0}^{2}\sqrt{2n_{\text{tot}}}x_{\mathrm{zp}}(\Omega_{0})\right.$
for the mechanical spring force. The criterion for the radiation pressure
force to exceed the mechanical spring force is established as
$\left.\sqrt{2}n_{\text{cav }}g_{0}>\sqrt{n_{\text{tot}}}\Omega_{0}\right.$.

In the context of our study on the membrane-in-the-middle optomechanical
device detailed in the main text, we find that shifting the resonance
frequency all the way to the zero requires only $9\times 10^{7}$ intracavity
photons. This photon number is lower than $4\times 10^{8}$ used in the
experiment in Ref. [36]. This comparison underscores the feasibility of our
approach in practical scenarios.

## References

## References

  * Riste  _et al._ [2013] D. Riste, M. Dukalski, C. A. Watson, G. De Lange, M. J. Tiggelman, Y. M. Blanter, K. W. Lehnert, R. N. Schouten, and L. DiCarlo, Nature 502, 350 (2013).
  * Cox  _et al._ [2016] K. C. Cox, G. P. Greve, J. M. Weiner, and J. K. Thompson, Physical Review Letters 116, 093602 (2016).
  * Sayrin  _et al._ [2011] C. Sayrin, I. Dotsenko, X. Zhou, B. Peaudecerf, T. Rybarczyk, S. Gleyzes, P. Rouchon, M. Mirrahimi, H. Amini, and M. Brune, Nature 477, 73 (2011).
  * Minev  _et al._ [2019] Z. K. Minev, S. O. Mundhada, S. Shankar, P. Reinhold, R. Gutiérrez-Jáuregui, R. J. Schoelkopf, M. Mirrahimi, H. J. Carmichael, and M. H. Devoret, Nature 570, 200 (2019).
  * Braginsky  _et al._ [1980] V. B. Braginsky, Y. I. Vorontsov, and K. S. Thorne, Science 209, 547 (1980).
  * Müller-Ebhardt  _et al._ [2009] H. Müller-Ebhardt, H. Rehbein, C. Li, Y. Mino, K. Somiya, R. Schnabel, K. Danzmann, and Y. Chen, Physical Review A 80, 043802 (2009).
  * Miao  _et al._ [2010] H. Miao, S. Danilishin, H. Müller-Ebhardt, and Y. Chen, New Journal of Physics 12, 083032 (2010).
  * Meng  _et al._ [2020] C. Meng, G. A. Brawley, J. S. Bennett, M. R. Vanner, and W. P. Bowen, Physical Review Letters 125, 043604 (2020).
  * de Lépinay  _et al._ [2021] L. M. de Lépinay, C. F. Ockeloen-Korppi, M. J. Woolley, and M. A. Sillanpää, Science 372, 625 (2021).
  * Kotler  _et al._ [2021] S. Kotler, G. A. Peterson, E. Shojaee, F. Lecocq, K. Cicak, A. Kwiatkowski, S. Geller, S. Glancy, E. Knill, and R. W. Simmonds, Science 372, 622 (2021).
  * Rossi  _et al._ [2018] M. Rossi, D. Mason, J. Chen, Y. Tsaturyan, and A. Schliesser, Nature 563, 53 (2018).
  * Guo  _et al._ [2019] J. Guo, R. Norte, and S. Gröblacher, Physical Review Letters 123, 223602 (2019).
  * Aspelmeyer  _et al._ [2014] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Reviews of Modern Physics 86, 1391 (2014).
  * Wilson  _et al._ [2015] D. J. Wilson, V. Sudhir, N. Piro, R. Schilling, A. Ghadimi, and T. J. Kippenberg, Nature 524, 325 (2015).
  * Meng  _et al._ [2022] C. Meng, G. A. Brawley, S. Khademi, E. M. Bridge, J. S. Bennett, and W. P. Bowen, Science Advances 8, eabm7585 (2022).
  * Bowen and Milburn [2015] W. P. Bowen and G. J. Milburn,  _Quantum Optomechanics_ (CRC Press, 2015).
  * Bennett  _et al._ [2016] J. S. Bennett, K. Khosla, L. S. Madsen, M. R. Vanner, H. Rubinsztein-Dunlop, and W. P. Bowen, New Journal of Physics 18, 053030 (2016).
  * Doherty  _et al._ [2012] A. C. Doherty, A. Szorkovszky, G. I. Harris, and W. P. Bowen, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 370, 5338 (2012).
  * Cohadon  _et al._ [1999] P. F. Cohadon, A. Heidmann, and M. Pinard, Physical Review Letters 83, 3174 (1999).
  * Arcizet  _et al._ [2006] O. Arcizet, P.-F. Cohadon, T. Briant, M. Pinard, and A. Heidmann, Nature 444, 71 (2006).
  * Kleckner and Bouwmeester [2006] D. Kleckner and D. Bouwmeester, Nature 444, 75 (2006).
  * Meng [2022] C. Meng,  _Optomechanical state preparation via strong measurement_ , PhD Thesis, The University of Queensland (2022).
  * Poggio  _et al._ [2007] M. Poggio, C. L. Degen, H. J. Mamin, and D. Rugar, Physical Review Letters 99, 017201 (2007).
  * Lee  _et al._ [2010] K. H. Lee, T. G. McRae, G. I. Harris, J. Knittel, and W. P. Bowen, Physical Review Letters 104, 123604 (2010).
  * Harris  _et al._ [2013] G. I. Harris, D. L. McAuslan, T. M. Stace, A. C. Doherty, and W. P. Bowen, Physical Review Letters 111, 103603 (2013).
  * Harris  _et al._ [2012] G. I. Harris, U. L. Andersen, J. Knittel, and W. P. Bowen, Physical Review A 85, 061802 (2012).
  * Whittle  _et al._ [2021] C. Whittle, E. D. Hall, S. Dwyer, N. Mavalvala, V. Sudhir, R. Abbott, A. Ananyeva, C. Austin, L. Barsotti, J. Betzwieser, C. D. Blair, A. F. Brooks, D. D. Brown, A. Buikema, C. Cahillane, J. C. Driggers, A. Effler, A. Fernandez-Galiana, P. Fritschel, V. V. Frolov, T. Hardwick, M. Kasprzack, K. Kawabe, N. Kijbunchoo, J. S. Kissel, G. L. Mansell, F. Matichard, L. McCuller, T. McRae, A. Mullavey, A. Pele, R. M. S. Schofield, D. Sigg, M. Tse, G. Vajente, D. C. Vander-Hyde, H. Yu, H. Yu, C. Adams, R. X. Adhikari, S. Appert, K. Arai, J. S. Areeda, Y. Asali, S. M. Aston, A. M. Baer, M. Ball, S. W. Ballmer, S. Banagiri, D. Barker, J. Bartlett, B. K. Berger, D. Bhattacharjee, G. Billingsley, S. Biscans, R. M. Blair, N. Bode, P. Booker, R. Bork, A. Bramley, K. C. Cannon, X. Chen, A. A. Ciobanu, F. Clara, C. M. Compton, S. J. Cooper, K. R. Corley, S. T. Countryman, P. B. Covas, D. C. Coyne, L. E. H. Datrier, D. Davis, C. Di Fronzo, K. L. Dooley, P. Dupej, T. Etzel, M. Evans, T. M. Evans, J. Feicht, P. Fulda, M. Fyffe, J. A. Giaime, K. D. Giardina, P. Godwin, E. Goetz, S. Gras, C. Gray, R. Gray, A. C. Green, E. K. Gustafson, R. Gustafson, J. Hanks, J. Hanson, R. K. Hasskew, M. C. Heintze, A. F. Helmling-Cornell, N. A. Holland, J. D. Jones, S. Kandhasamy, S. Karki, P. J. King, R. Kumar, M. Landry, B. B. Lane, B. Lantz, M. Laxen, Y. K. Lecoeuche, J. Leviton, J. Liu, M. Lormand, A. P. Lundgren, R. Macas, M. MacInnis, D. M. Macleod, S. Márka, Z. Márka, D. V. Martynov, K. Mason, T. J. Massinger, R. McCarthy, D. E. McClelland, S. McCormick, J. McIver, G. Mendell, K. Merfeld, E. L. Merilh, F. Meylahn, T. Mistry, R. Mittleman, G. Moreno, C. M. Mow-Lowry, S. Mozzon, T. J. N. Nelson, P. Nguyen, L. K. Nuttall, J. Oberling, R. J. Oram, C. Osthelder, D. J. Ottaway, H. Overmier, J. R. Palamos, W. Parker, E. Payne, R. Penhorwood, C. J. Perez, M. Pirello, H. Radkins, K. E. Ramirez, J. W. Richardson, K. Riles, N. A. Robertson, J. G. Rollins, C. L. Romel, J. H. Romie, M. P. Ross, K. Ryan, T. Sadecki, E. J. Sanchez, L. E. Sanchez, T. R. Saravanan, R. L. Savage, D. Schaetz, R. Schnabel, E. Schwartz, D. Sellers, T. Shaffer, B. J. J. Slagmolen, J. R. Smith, S. Soni, B. Sorazu, A. P. Spencer, K. A. Strain, L. Sun, M. J. Szczepańczyk, M. Thomas, P. Thomas, K. A. Thorne, K. Toland, C. I. Torrie, G. Traylor, A. L. Urban, G. Valdes, P. J. Veitch, K. Venkateswara, G. Venugopalan, A. D. Viets, T. Vo, C. Vorvick, M. Wade, R. L. Ward, J. Warner, B. Weaver, R. Weiss, B. Willke, C. C. Wipf, L. Xiao, H. Yamamoto, L. Zhang, M. E. Zucker, and J. Zweizig, Science 372, 1333 (2021).
  * Note [1] $\chi_{\mathrm{fb}}^{-1}(\omega)=-Km(1-i\omega/\Gamma_{\mathrm{fb}})e^{i\omega\tau}\approx-Km(1+\omega^{2}\tau/\Gamma_{\mathrm{fb}}+i\omega(\tau-1/\Gamma_{\mathrm{fb}}))$.
  * Buonanno and Chen [2002] A. Buonanno and Y. Chen, Physical Review D 65, 042001 (2002).
  * Chen  _et al._ [2023] J. Chen, B. B. Lane, S. Direkci, D. Ganapathy, X. Yin, N. Mavalvala, Y. Chen, and V. Sudhir, Causal State Estimation and the Heisenberg Uncertainty Principle (2023), arXiv:2304.14476 [quant-ph].
  * Note [100] Supplemental material.
  * Wiener [1964] N. Wiener,  _Extrapolation, Interpolation, and Smoothing of Stationary Time Series_ (The MIT Press, 1964).
  * Kohen  _et al._ [1997] D. Kohen, C. C. Marston, and D. J. Tannor, The Journal of Chemical Physics 107, 5236 (1997).
  * Note [2] Inside the backaction dominated regime, the squeezing criterion $C>\left(n_{\mathrm{tot}}Q^{2}/4\right)^{1/3}$ can be simplified to $\left.C>Q/2\right.$.
  * Khosla  _et al._ [2017] K. E. Khosla, G. A. Brawley, M. R. Vanner, and W. P. Bowen, Optica 4, 1382 (2017).
  * Purdy  _et al._ [2013] T. P. Purdy, R. W. Peterson, and C. A. Regal, Science 339, 801 (2013).
  * Wilson [2012] D. J. Wilson,  _Cavity Optomechanics with High Stress Silicon Nitride Films_ , phd, California Institute of Technology (2012).
  * Clerk  _et al._ [2010] A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, Reviews of Modern Physics 82, 1155 (2010).
  * Note [3] The protocol does not accommodate the weak measurement regime, as the thermal motion cannot be resolved from the measurement, thereby rendering feedback application unfeasible. Furthermore, it is worth noting that, through feedback control, shifting $\Omega$ down nearly to DC only requires $\Omega_{\mathrm{meas}}\gg\sqrt{\Omega_{0}\Gamma}$ while ground state cooling requires $\Omega_{\mathrm{meas}}\gg 2n_{\mathrm{th}}n_{\mathrm{tot}}\sqrt{\Omega_{0}\Gamma}$, i.e., $\eta C\gg n_{\mathrm{th}}$. Therefore, feedback softening does not have as significant a technical hurdle as feedback cooling.