Is there any way how to compute convolution matrix for Nearest Neighbor (bicubic, bilinear) image scaling (upscale/downscale)?
This can't be done as a convolution because some of the operations you mentioned are not linear. Additionally, convolution by definition (multiplication by a toplitz matrix) shouldn't ever change the input and output size.
Notice the high order terms.
As others have mentioned, convolution operation cannot change the size of the image. However, interpolation is a convolution operation. You just need to upsample by zero padding the image before performing the convolution (filtering) operation.
For example, say you want to interpolate a 1D signal times 2. The original signal is given by $a_0, a_1, .. , a_n$. You first zero pad it to get $a_0, 0, a_1, 0, ..$, then you filter the result. If the filter taps are
0.5, 1, 0.5, you will get simple linear interpolation. Similar scheme can be used in 2D images for biliniear interpolation, spline interpolation, etc.