Here is an attempt, tell me what I misunderstood.
You say that $F\times F^{H}$ (which is of dimension $m\times m$ ) is unitary, which implies it is invertible. Therefore it means that $F$ is of full rank, i.e $rank(F) = m$.
From that we know from the rank nullity theorem that $ker(F)$ will be of dimension $N-m$. This means that no matter how hard you try, you won't find more than $N-m$ independent vectors $d_i$ such that $F \times d_i = 0$. This seriously compromise the existence of a G of size $N\times N$. Instead it will be of size $N\times(N-m)$.
Now how to find such a G matrix? For what I understand of singular value decomposition it exists no matter what are the dimensions of your matrix. In particular it exists if $m<N$. So we can write:
$$F = U\Sigma V^{H}\tag{1}\label{eq1}$$
with $U$ of dimension $m\times m$ unitary, $V$ of dimension $N\times N$ and also unitary. We know from the fullness of $F$'s rank that
$$\Sigma = \begin{bmatrix} diag(\lambda_1, ..., \lambda_m) && 0_{m\times (N-m)}\end{bmatrix}\tag{2}\label{eq2}$$
with $\lambda_1$, ..., $\lambda_m$ all nonzero.
Now if you multiply $\Sigma$ by :
$$e_{m+p} = \begin{bmatrix}0\\0\\...\\0\\ 1\\0\\...\\0 \end{bmatrix}\tag{3}\label{eq3}$$
with the $1$ at the n+p-th position and $1\leq p\leq N-m$ we get $0$. From that with
$$G = V\times \begin{bmatrix}e_{m+1}, e_{m+2}, ..., e_N\end{bmatrix}\tag{4}\label{eq4}$$
you get
$$F\times G = U\Sigma V^{H} \times V\times\begin{bmatrix}e_{m+1}, e_{m+2}, ..., e_N\end{bmatrix} = U \begin{bmatrix}\Sigma e_{m+1},\Sigma e_{m+2}, ..., \Sigma e_N\end{bmatrix} = 0 \tag{5}\label{eq5}$$
One can easily verify :
$$(Ve_{n+i})^{H}\times (Ve_{n+j}) = e_{n+i}^Te_{n+j} = \delta_{ij} \tag{6}\label{eq6}$$
So that the $Ve_{n+k}$ with $1\leq k\leq N-m$ form an orthonormal basis of $ker(F)$ which is the closest thing you can get from the G you specified.
I also found something else :
Since you said $F\times F^{H}$ is unitary, using the SVD decomposition from \ref{eq1} the unitary nature of $F\times F^{H}$ we get :
$$F\times F^{H}\times F\times F^{H} = U\Sigma V^{H}V\Sigma^{H} U^{H}U\Sigma V^{H}V\Sigma^{H} U^{H} = U (\Sigma\Sigma^{H})^{2} U^{H} = I_{m\times m}\tag{7}\label{eq7}$$
This yields $$(\Sigma\Sigma^{H})^{2}= I_{m\times m}\tag{8}\label{eq8}$$
from \ref{eq2} we easily get $$\Sigma\Sigma^{H} = diag(\vert\lambda_1\vert^{2}, \vert\lambda_2\vert^{2}, ..., \vert\lambda_m\vert^{2})\tag{9}\label{eq9}$$ Therefore \ref{eq7} translates into :
$$diag(\vert \lambda_1\vert^{4}, \vert \lambda_2\vert^{4}, ..., \vert \lambda_m\vert^{4}) = I_{m\times m}\tag{10}\label{eq10}$$
which finally means $$\vert\lambda_i\vert =1\quad\forall 1\leq i\leq m\tag{11}\label{eq11}$$
I found it interesting.
If this does not answer your question please tell us what you mean by "at least the result is minimum" (which result? minimum in what regards?). Also could you detail what you mean by "can be found directly or it must to be optimized"? Anyway I hope it helped you one way or another.
Edit :
You seem to be interested in convergence properties so I went back looking at how SVD works (I had forgotten). In a naive algorithm (like if I had to implement it from scrathes myself today) I would have to compute $A\times A^{H}$, then the eigen elements of this matrix (using the so-called QR algorithm, which is a converging and to my knowledge pretty stable algorithm).
The eigenvectors would give me the U matrix, Then I'd multiply $A^{H}$ with some of those eigenvectors (those associated with nonzero eigenvalues) and complete this into a orthonormal basis of $\mathbb{R}^m$ and this would give me V. and the last $m-r$ (where $r = rank(A)$) vectors of V would give G.
So the most numerically perilous parts here seem to be the computation of $A\times A^{H}$ and then the computation of the $A^{H} U_i$ (with $U_i$ the ad hoc columns of $U$).
That being said here is a code made in python checking the computation of G using numpy's SVD and the errors made computing $A\times G$.
I simplified a bit, multiplying by the e matrix is equivalent to selecting its last (N-m) columns
import numpy as np
n = 100
m = 1000
A = np.random.rand(n, m)
svd_decomposition = np.linalg.svd(A)
U = svd_decomposition[0]
V = svd_decomposition[2].T
r = np.linalg.matrix_rank(A)
G = V[:, r:]
res = np.abs(A@G).sum()
print(res)
with the columns and row numbers I set I get something in the 10-12 realm. The error seem to grow with the number of rows but is pretty insensitive to the number of columns.
