  We develop a methodology for the computation of the $K\to \ell\nu_\ell
\ell'^+ \ell'^-$ decay width using lattice QCD and present an exploratory study
here. We use a scalar function method to account for the momentum dependence of
the decay amplitude and adopt the infinite volume reconstruction (IVR) method
to reduce the systematic errors such as the temporal truncation effects and the
finite-volume effects. We then perform a four-body phase-space integral to
obtain the decay width. The only remaining technical problem is the possible
power-law finite-volume effects associated with the process of $K\to\pi\pi
\ell\nu_\ell\to \ell\nu_\ell \ell'^+ \ell'^-$, where the intermediate state
involves multiple hadrons. In this work, we use a gauge ensemble of twisted
mass fermion with a pion mass $m_\pi=352$ MeV and a nearly-physical kaon mass.
At this kinematics, the $\pi\pi$ in the intermediate state cannot be on shell
simultaneously as $2m_\pi>m_K$ and the finite-volume effects associate with
$\pi\pi$ state are exponentially suppressed. Using the developed methods
mentioned above, we calculate the branching ratios for four channels of $K\to
\ell\nu_\ell\ell'^+ \ell'^-$, and obtain the results close to the experimental
measurements and ChPT predictions. Our work demonstrates the capability of
lattice QCD to improve Standard Model prediction in $K\to \ell\nu_\ell \ell'^+
\ell'^-$ decay width.
