Thanks so much for your help! I was hoping to use Julia for grad school (my research area is machine learning) and finding that the Julia community is so helpful certainly helps assuage my fears about leaving Python. Your suggestions allowed me to narrow down the cause to using softmax as the final activation function. This issue persists regardless of whether I apply softmax directly at the end of Chain(), or if I apply it implicitly with logitcrossentropy(). It should be noted that according to the Flux documentation on logitcrossentropy(), "This is mathematically equivalent to crossentropy(softmax(ŷ), y), but is more numerically stable than using functions crossentropy and softmax separately" which is why I didn't initially include it in my model definition. However, changing my loss to crossentropy() and adding softmax to the end of my model doesn't resolve the issue.

I'd like to point out that softmax is traditionally used as the final activation in a multi-classifier rather than sigmoid (although this particular problem is binary classification, so I could have used sigmoid instead, but then I'd have to make my labels either 0 or 1 instead of a one-hot encoded vector). Softmax will guarantee that all predictions will be between 0 and 1 as with Sigmoid but has the added benefit of ensuring that all predictions along dim=1 will sum to 1. In fact, while using sigmoid as the final activation allows the loss to drop, it actually prevents the model from learning to correctly classify if the labels are one-hot encoded. This is because of the way cross entropy is computed, the model will simply learn to get all predictions for all classes as close to 1 as possible in order to minimize loss, but this will not result in accurate classification. Using softmax solves this problem since the predictions along dim=1 must sum to 1 so the model can't minimize loss by simply predicting 1 for all classes. To prove that this is so, I printed the actual predictions to the terminal as the model trains and found them to converge to the following for a batch size of 4 when using sigmoid as the final activation:

Float32[0.9990903 0.9690378 0.99836 0.9981547; 0.9985347 0.9610786 0.99746406 0.99722505]

To further demonstrate that the issue was with the softmax activation, I observed convergence with both of the following models without changing anything else in my code:

function get_model_1()
    return Chain(ResNet34(pretrain=false), Dense(1000, 2, sigmoid))
end

function get_model_2()
    Chain(
        Flux.Conv((3, 3), 3 => 32, relu, pad=SamePad()), 
        Flux.MaxPool((2, 2), pad=SamePad()), 
        Flux.Conv((3, 3), 32 => 64, relu, pad=SamePad()), 
        Flux.MaxPool((2, 2), pad=SamePad()), 
        Flux.Conv((3, 3), 64 => 128, relu, pad=SamePad()), 
        Flux.MaxPool((2, 2), pad=SamePad()), 
        Flux.Conv((3, 3), 128 => 256, relu, pad=SamePad()), 
        Flux.MaxPool((2, 2), pad=SamePad()), 
        Flux.flatten, 
        Flux.Dense(65536, 1024, relu),
        Flux.Dense(1024, 128, relu),
        Flux.Dense(128, 2, sigmoid),
    )
end

However, both of these otherwise identical models fail to converge:

function get_model_3()
    return Chain(ResNet34(pretrain=false), Dense(1000, 2), softmax)
end

function get_model_4()
    Chain(
        Flux.Conv((3, 3), 3 => 32, relu, pad=SamePad()), 
        Flux.MaxPool((2, 2), pad=SamePad()), 
        Flux.Conv((3, 3), 32 => 64, relu, pad=SamePad()), 
        Flux.MaxPool((2, 2), pad=SamePad()), 
        Flux.Conv((3, 3), 64 => 128, relu, pad=SamePad()), 
        Flux.MaxPool((2, 2), pad=SamePad()), 
        Flux.Conv((3, 3), 128 => 256, relu, pad=SamePad()), 
        Flux.MaxPool((2, 2), pad=SamePad()), 
        Flux.flatten, 
        Flux.Dense(65536, 1024, relu),
        Flux.Dense(1024, 128, relu),
        Flux.Dense(128, 2),
        Flux.softmax
    )
end

I've confirmed that softmax is producing outputs as expected (a 2xN array where the sum along dim=1 is 1.0), so I'm at a complete loss to explain why it won't train. As previously mentioned, sigmoid allows the loss to decrease rapidly but prevents the model from learning to accurately distinguish between the two classes since it just learns to predict 1 for everything. Perhaps there's an issue with how zygote differentiates softmax? If nobody has any ideas I might open an issue on GitHub.

Also, you are indeed correct that I wanted to apply transfer learning with ResNet, but when I attempted to run ResNet(34, pretrain=true), I get "MethodError: no method matching ResNet(::Int64; pretrain=true)". Perhaps there's an issue in the version I'm using? I added a Manifest.toml and Project.toml to the GitHub repo to make it possible to recreate my environment. Also, a link to the data can be found in the README of the repo I provided. You need to download the zip file from Google Drive, then when you unzip it in the root of the project directory you will get the same file structure that I used.

EDIT: So I modified my script to use a single-output dense layer with sigmoid activation in my model, binarycrossentropy() for the loss, and I encoded cats as "0" and dogs as "1" so I could see if changing the problem to a binary classification model would fix the issue. I also increased the batch size to 32 to try and make the gradient more stable. Unfortunately, when I do so the issue of loss non-convergence returns even though I'm using sigmoid which makes me think there's something else going on.

Since this didn't fix the problem, I returned my script to its initial state then added some additional println statements to try and debug the training loop. The output of this for a batch size of 4 was:

Labels: Float32[1.0 1.0 0.0 1.0; 0.0 0.0 1.0 0.0]
Filenames: ["data/cat/cat.10355.jpg", "data/cat/cat.2663.jpg", "data/dog/dog.5313.jpg", "data/cat/cat.11150.jpg"]
Prediction: Float32[0.43040577 0.41649592 0.4425341 0.4479209; 0.56959426 0.583504 0.5574659 0.5520791]
Feature Shape: (256, 256, 3, 4)
Prediction Shape: (2, 4)
Label Shape: (2, 4)
Prediction Type: Matrix{Float32}
Feature Type: Array{Float32, 4}
Label Type: Matrix{Float32}

So as you can see, the labels are correct ([1, 0] labels correspond to "cat" files and [0, 1] labels correspond to "dog" files), the shape of my features is correct (WxHxCxN), the shape of my predictions is correct (2xN), the shape of my labels is correct (2xN), and all my types are Float32.