You are right in saying that $I$ allows one to relate angular velocity and angular momentum in a linear way. It is just not as simple as the momentum and velocity case. An intuition for why things get complicated is that $L = r \times p $ involves a cross product which makes it very sensitive to the choice of a specific set of orthonormal bases(with fixed origin). While $p=mv$ involves a scalar mass that is independent of your choice of coordinates.
To explain inertia tensor, I guess we could start with simpler cases where sufficient symmetry is present (for example a sphere in 3D or a circular pancake in 2d), $L = I(\omega)$ reduces to $L = I\omega$ where $L$ and $\omega$ are vectors and $I$ is just a scalar. An intuition for this reduction is that symmetry makes $I$ resemble $m$ more. As I mentioned earlier, mass is always independent of coordinate choice. but $I$ is only independent of the coordinates that PRESERVE symmetry. Therefore, spheres and circular pancakes are pretty easy to deal with and no inertia tensor is necessary.
But for a general, extended, rigid body in 3d, the lack of symmetry breaks the simple linear relationship. Suppose you have an orthonormal basis, the origin of which is the corner of a cube and the axes line up with the edges of the cube. Basically when the cube is rotated around the z-axis, all the parts of the cube are also instantaneously rotating in the other directions (if you draw a diagram, it would be clear). Therefore $\omega_z$ affects $L_y$ and $L_x$. I don't see an intuitive explanation of the quantitative details.. But this simple cube example shows that $L_x, L_y, L_z$ must each be a linear combination of $\omega_x, \omega_y, \omega_z$. And the mathematical expression that quantifies this must be a matrix.
*A mathematical sidetrack: this matrix itself is not a tensor, but rather a REPRESENTATION of a tensor that maps angular velocity vectors to angular momentum DUAL vectors. In abstract index notation, $L_\alpha = I_{\alpha\beta}\omega^\beta $ You will see a lot of similar notations in E&M, Relativity etc.
I found in some game physics engines the equation is expressed in inertia frame rather than body frame.
$$\tau=I\dot{\omega}+\omega \times I \omega$$
In that situation the inertia tensor $I$ is not constant but varying with the body frame. Therefore it is updated with the orientation of rigid body. Moment $\tau$ and angular velocity $\omega$ are expressed in world frame(inertia frame). Is the formula above established?
tl;dr: Yes, this is valid. Getting there is nontrivial.
The $\boldsymbol \omega_B \times (\mathrm I_B \, \boldsymbol \omega_B)$ term in $\boldsymbol \tau_B = \mathrm I_B \, \dot{\boldsymbol \omega}_B + \boldsymbol \omega_B \times (\mathrm I_B \, \boldsymbol \omega_B)$ results from working in a rotating frame. It is a fictitious torque, very much akin to the fictitious forces that arise with Newton's second law when computed in a rotating frame. Just as those fictitious forces vanish in an inertial frame, so does the fictitious torque. The inertial frame rotational equation of motion is
$$\boldsymbol\tau_I = \frac d{dt} (\mathrm I_I\, \boldsymbol \omega_I) =
\mathrm I_I \, \dot {\boldsymbol \omega}_I
+ \dot {\mathrm I}_I \,\boldsymbol \omega_I$$
This raises the question, what is the time derivative of the inertia tensor from the perspective of an inertial frame? To determine this, the relationship between the inertia tensor as expressed in the body frame and as expressed inertial frame is needed. This is
$$\mathrm I_I = \mathrm T_{R\to I}\,\mathrm I_B\,\mathrm T_{R\to I}^{\,T}$$
where $\mathrm I_I$ is the inertia tensor expressed in the inertial frame, $\mathrm T_{R\to I}$ is the transformation matrix from the body frame to the inertial frame, $\mathrm I_B$ the inertia tensor expressed in the body frame, and $\mathrm T_{R\to I}^{\,T}$ is the transpose of (and also the inverse of) $\mathrm T_{R\to I}$. Taking the derivative with respect to time yields
$$\dot {\mathrm I}_I =
\dot {\mathrm T}_{R\to I} \, \mathrm I_B\,\mathrm T_{R\to I}^{\,T}
+ \mathrm T_{R\to I} \, \mathrm I_B \, \dot {\mathrm T}_{R\to I}^{\,T}$$
Without derivation, the time derivative of the transformation matrix is
$$\dot {\mathrm T}_{R\to I} = \mathrm T_{R\to I} \operatorname{Sk}(\boldsymbol \omega_B) = \operatorname{Sk}(\boldsymbol \omega_I)\mathrm T_{R\to I}$$
where $\operatorname{Sk}(a)$ is the skew symmetric cross product matrix generated by the vector $\boldsymbol a$:
$$\operatorname{Sk}(\boldsymbol a) \equiv \begin{pmatrix}
\phantom{-}0 & -a_z & \phantom{-}a_y \\
\phantom{-}a_z & \phantom{-}0 & -a_x \\
-a_y & \phantom{-}a_x & \phantom{-}0
\end{pmatrix}$$
Note that $\operatorname{Sk}(\boldsymbol a) \, \boldsymbol b = \boldsymbol a \times \boldsymbol b$. With the above, the time derivative of the inertia tensor as expressed in the inertial frame is
$$\dot {\mathrm I}_I =
\operatorname{Sk}(\boldsymbol \omega_I) \mathrm T_{R\to I} \, \mathrm I_B\,\mathrm T_{R\to I}^{\,T}
+ \mathrm T_{R\to I} \, \mathrm I_B \, \mathrm T_{R\to I}^{\,T} \operatorname{Sk}^T(\boldsymbol \omega_I)
= \operatorname{Sk}(\boldsymbol \omega_I)\,\mathrm I_I - \mathrm I_I \, \operatorname{Sk}(\boldsymbol \omega_I)$$
This can be nonzero due to the non-commutative nature of matrix multiplication. We're not interested so much in the form of $\dot {\mathrm I}_I$ in general so much as the product $\dot {\mathrm I}_I\,\boldsymbol \omega_I$. This is
$$\dot {\mathrm I}_I\,\boldsymbol \omega_I =
\operatorname{Sk}(\boldsymbol \omega_I)\,\mathrm I_I \,\boldsymbol \omega_I
- \mathrm I_I \, \operatorname{Sk}(\boldsymbol \omega_I) \,\boldsymbol \omega_I
= \boldsymbol \omega_I \times (\mathrm I_I \,\boldsymbol \omega_I)
- \mathrm I_I \, (\boldsymbol \omega_I \times \boldsymbol \omega_I)
= \boldsymbol \omega_I \times (\mathrm I_I \,\boldsymbol \omega_I)$$
and thus
$$\boldsymbol\tau_I =
\mathrm I_I \, \dot {\boldsymbol \omega}_I
+ \boldsymbol \omega_I \times (\mathrm I_I \,\boldsymbol \omega_I)$$
Best Answer
Translation
\begin{align*} &\mathbf{v}_I=\mathbf{S}\,\mathbf{v}_B\quad\Rightarrow\, \mathbf{\dot{v}}_I=\mathbf{\dot S}\,\mathbf{v}_B+\mathbf{S}\,\mathbf{\dot{v}}_B\\ &\text{with}\quad \mathbf{\dot S}=\mathbf{S}\,\mathbf{\omega}_B^\times\\ &\Rightarrow \end{align*} \begin{align*} &\mathbf{\dot{v}}_I=\mathbf{S}\,\left(\mathbf{\omega}_B\times \,\mathbf{v}_B\right)+\mathbf{S}\,\mathbf{\dot{v}}_B\\ \end{align*} Newton equation \begin{align*} &m\,\mathbf{\dot{v}}_I=\mathbf{F}_I\quad \text{or}\quad, m\,\underbrace{\mathbf{S}^T\,\mathbf{\dot{v}}_I}_{(\mathbf{\dot{v}}_I)_B}=\mathbf{S}^T\,\mathbf{F}_I=\mathbf{F}_B \end{align*} you obtain \begin{align*} &\mathbf{S}^T\,\mathbf{\dot{v}}_I=\left(\mathbf{\omega}_B\times \,\mathbf{v}_B\right)+\mathbf{\dot{v}}_B=\frac{\mathbf{F}_B}{m}\\ &\boxed{\,\mathbf{\dot{v}}_B=\frac{\mathbf{F}_B}{m}-\left(\mathbf{\omega}_B\times \,\mathbf{v}_B\right)\,} \end{align*}
Rotation
Euler equation in B_system
\begin{align*} &\boxed{\,I_B\,\mathbf{\dot{\omega}}_B+\mathbf{\omega}_B\,\times \left(I_B\,\mathbf{\omega}_B\right)=\mathbf\tau_B\,} \end{align*}
Notice that with those two equations of motion, you don't get the position and the orientation (angles) of the rigid body , you need additional equations
Edit
Rotation Matrix $~\mathbf{S}~$
\begin{align*} &\mathbf{S}=\mathbf{S}_z(\psi)\,\mathbf{S}_x(\varphi)\,\mathbf{S}_y(\vartheta)\\ &\mathbf{S}= \left[ \begin {array}{ccc} \cos \left( \psi \right) &-\sin \left( \psi \right) &0\\ \sin \left( \psi \right) &\cos \left( \psi \right) &0\\ 0&0&1\end {array} \right] \, \left[ \begin {array}{ccc} 1&0&0\\ 0&\cos \left( \varphi \right) &-\sin \left( \varphi \right) \\ 0 &\sin \left( \varphi \right) &\cos \left( \varphi \right) \end {array} \right] \, \left[ \begin {array}{ccc} \cos \left( \vartheta \right) &0&\sin \left( \vartheta \right) \\ 0&1&0 \\ -\sin \left( \vartheta \right) &0&\cos \left( \vartheta \right) \end {array} \right] \\ &\text{with}\quad \mathbf{\dot S}=\mathbf{S}\,\mathbf{\omega}^\times\\ &\Rightarrow\\ &\begin{bmatrix} \omega_x \\ \omega_y \\ \omega_z \\ \end{bmatrix} =\underbrace{ \left[ \begin {array}{ccc} \cos \left( \vartheta \right) &0&-\cos \left( \varphi \right) \sin \left( \vartheta \right) \\ 0&1&\sin \left( \varphi \right) \\ \sin \left( \vartheta \right) &0&\cos \left( \varphi \right) \cos \left( \vartheta \right) \end {array} \right] }_{\mathbf{J}_R} \,\begin{bmatrix} \dot{\varphi} \\ \dot{\vartheta} \\ \dot{\psi} \\ \end{bmatrix}\\ &\Rightarrow \end{align*} \begin{align*} & \boxed{\,\begin{bmatrix} \dot{\varphi} \\ \dot{\vartheta} \\ \dot{\psi} \\ \end{bmatrix}=\left[ \begin {array}{ccc} \cos \left( \vartheta \right) &0&\sin \left( \vartheta \right) \\ {\frac {\sin \left( \varphi \right) \sin \left( \vartheta \right) }{\cos \left( \varphi \right) }}&1&-{\frac {\sin \left( \varphi \right) \cos \left( \vartheta \right) }{\cos \left( \varphi \right) }} \\ -{\frac {\sin \left( \vartheta \right) }{\cos \left( \varphi \right) }}&0&{\frac {\cos \left( \vartheta \right) } {\cos \left( \varphi \right) }}\end {array} \right] \, \begin{bmatrix} \omega_x \\ \omega_y \\ \omega_z \\ \end{bmatrix}\,}\tag 3 \end{align*}
Singularity at $~\varphi=\pi/2~$
Inertial Position vector $~\mathbf R_I~$
\begin{align*} &\mathbf{v}_I=\mathbf{S}\,\mathbf{v}_B\\&\Rightarrow \end{align*} \begin{align*} & \boxed{\,\mathbf{\dot{R}}_I=\mathbf{v}_I\,}\tag 4 \end{align*}
all together you obtained 12 first order differential equations for a rigid body solution