Origin of Virtual Molecular Orbitals in Hartree-Fock equations

Let's look at scheme

Do I understand correctly, the origin of Virtual Orbitals is a solution of (for, exampre, restricted) HF equation eigenfuncions?

Yes, they result from the HF equations just like the occupied orbitals. The diagram should also be show them as present in the initial set of equations, unless the basis (used for the LCAO) was changed between the two states shown. The number of occupied orbitals plus the number of virtual orbitals is equal to the basis set size because they are the eigensolutions of $$\mathbf{F}\mathbf{C} = \mathbf{\epsilon}\mathbf{S}\mathbf{C},$$ where $$\mathbf{F}$$ is the Fock matrix, $$\mathbf{C}$$ is the matrix of orbital coefficients, $$\mathbf{\epsilon}$$ are the eigenvalues which are the orbital energies, and $$\mathbf{S}$$ is the atomic orbital overlap matrix (because the basis set used is typically not orthogonal). The dimension of all matrices involved is determined by the size of the basis set, which will be finite because otherwise the computational effort becomes infinite.